diff --git a/Project.toml b/Project.toml index 2af65cc..1eee064 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "PowerSystemsMaps" uuid = "e146f72c-d4f8-45d3-b3ca-12a16cb68a38" authors = ["cbarrows "] -version = "0.2.0" +version = "0.2.1" [deps] Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" diff --git a/src/PowerSystemsMaps.jl b/src/PowerSystemsMaps.jl index a9da121..78ee280 100644 --- a/src/PowerSystemsMaps.jl +++ b/src/PowerSystemsMaps.jl @@ -15,8 +15,8 @@ export plot_net export plot_net! export plot_components! export make_graph +export plot_map include("plot_network.jl") -include("sfdp_fixed.jl") end # module PowerSystemsMaps diff --git a/src/plot_network.jl b/src/plot_network.jl index 56d8121..2729eb2 100644 --- a/src/plot_network.jl +++ b/src/plot_network.jl @@ -82,12 +82,10 @@ function make_graph(sys::PowerSystems.System; kwargs...) :name => get_name(b), :number => get_number(b), :area => get_name(get_area(b)), - :fixed => false, ) if has_supplemental_attributes(GeographicInfo, b) (lon, lat) = get_lonlat(b) data[:initial_position] = PT(lon, lat) - data[:fixed] = true end g[get_name(b)] = data end @@ -105,24 +103,22 @@ function make_graph(sys::PowerSystems.System; kwargs...) # layout nodes a = adjacency_matrix(g) # generates a sparse adjacency matrix - fixed = get_prop(g, :fixed) - sort_ids = vcat(findall(fixed), findall(.!fixed)) - orig_ids = sortperm(sort_ids) - sorted_a = a[sort_ids, sort_ids] - - @info "calculating node locations with SFDP_Fixed" - K = get(kwargs, :K, 0.1) - ip = get_prop(g, :initial_position) - setdiff!(ip, ip[findall(isnothing, ip)]) - network = sfdp_fixed( - sorted_a; - tol = 1.0, - C = 0.0002, - K = K, - iterations = 100, - fixed = true, - initialpos = ip, - )[orig_ids] # generate 2D layout and sort back to order of a + ip = Dict(zip(1:nv(g), get_prop(g, :initial_position))) + filter!(p -> !isnothing(last(p)), ip) + + if length(ip) != nv(g) + @info "calculating node locations with SFDP" + network = NetworkLayout.sfdp( + a; + tol = get(kwargs, :tol, 1.0), + C = get(kwargs, :C, 0.0002), + K = get(kwargs, :K, 0.1), + iterations = get(kwargs, :iterations, 100), + pin = ip, + ) + else + network = ip + end set_prop!(g, :x, first.(network)) set_prop!(g, :y, last.(network)) @@ -256,12 +252,64 @@ function plot_net!(p::Plots.Plot, g; kwargs...) return p end +""" +map network on top of basemap + +Accepted kwargs for network plot: + - `lines::Bool` : show lines + - `linealpha::Union{Float, Vector{Float}}` : line transparency + - `nodesize::Union{Float, Vector{Float}}` : node size + - `nodehover::Union{String, Vector{String}}` : node hover text + - `linewidth::Float` : width of lines + - `linecolor::String` : line color + - `buffer::Float` + - `size::Tuple` : figure size + - `xlim::Tuple` : crop x axis + - `ylim::Tuple` : crop y axis + - `showleged::Bool` : show legend + - `nodecolor::String` : node color + - `nodealpha::Float` : node transparency + - `legend_font_color::String` : legend font color + +Accepted kwargs for map: + - `map_lines::Bool` : show lines + - `map_linealpha::Union{Float, Vector{Float}}` : line transparency + - `map_nodesize::Union{Float, Vector{Float}}` : node size + - `map_nodehover::Union{String, Vector{String}}` : node hover text + - `map_linewidth::Float` : width of lines + - `map_linecolor::String` : line color + - `map_buffer::Float` + - `map_size::Tuple` : figure size + - `map_xlim::Tuple` : crop x axis + - `map_ylim::Tuple` : crop y axis + - `map_showleged::Bool` : show legend + - `map_nodecolor::String` : node color + - `map_nodealpha::Float` : node transparency + - `map_legend_font_color::String` : legend font color +""" +function plot_map(sys::System, shapefile::AbstractString; kwargs...) + g = make_graph(sys; kwargs...) + shp = Shapefile.shapes(Shapefile.Table(shapefile)) + shp = lonlat_to_webmercator(shp) #adjust coordinates + + map_kwargs = filter(x -> startswith(string(first(x)), "map_"), kwargs) + kwargs = setdiff(kwargs, map_kwargs) + map_kwargs = [Symbol(replace(string(k), "map_" => "")) => v for (k, v) in map_kwargs] + + # plot a map from shapefile + p = plot(shp; map_kwargs...) + + # plot the network on the map + p = plot_net!(p, g; kwargs...) + return p +end + # borrowed from function lonlat_to_webmercator(xLon, yLat) # Check coordinates are in range - abs(xLon) <= 180 || throw("Maximum longitude is 180.") - abs(yLat) < 85.051129 || throw( + abs(xLon) > 180 && throw("Maximum longitude is 180.") + abs(yLat) >= 85.051129 && throw( "Web Mercator maximum lattitude is 85.051129. This is the lattitude at which the full map becomes a square.", ) diff --git a/src/sfdp_fixed.jl b/src/sfdp_fixed.jl deleted file mode 100644 index 812259a..0000000 --- a/src/sfdp_fixed.jl +++ /dev/null @@ -1,188 +0,0 @@ -#export SFDP_Fixed, sfdp_fixed - -#extras -const NL = NetworkLayout -import NetworkLayout: norm -sfdp_fixed(g; kwargs...) = NL.layout(SFDP_Fixed(; kwargs...), g) - -""" - SFDP_Fixed(; kwargs...)(adj_matrix) - sfdp_fixed(adj_matrix; kwargs...) - -Using the Spring-Electric [model suggested by Yifan Hu](http://yifanhu.net/PUB/graph_draw_small.pdf). -Forces are calculated as: - - f_attr(i,j) = ‖xi - xj‖ ² / K , i<->j - f_repln(i,j) = -CK² / ‖xi - xj‖ , i!=j - -Takes adjacency matrix representation of a network and returns coordinates of -the nodes. - -## Keyword Arguments -- `dim=2`, `Ptype=Float64`: Determines dimension and output type `NL.Point{dim,Ptype}`. -- `tol=1.0`: Stop if position changes of last step `Δp <= tol*K` for all nodes -- `C=0.2`, `K=1.0`: Parameters to tweak forces. -- `iterations=100`: maximum number of iterations -- `initialpos=NL.Point{dim,Ptype}[]` - - Provide list of initial positions. If length does not match Network size the initial - positions will be truncated or filled up with random values between [-1,1] in every coordinate. - -- `seed=1`: Seed for random initial positions. -- `fixed = false` : Anchors initial positions. -""" -struct SFDP_Fixed{Dim, Ptype, T <: AbstractFloat} <: IterativeLayout{Dim, Ptype} - tol::T - C::T - K::T - iterations::Int - initialpos::Vector{NL.Point{Dim, Ptype}} - seed::UInt - fixed::Bool -end - -# TODO: check SFDP_Fixed default parameters -function SFDP_Fixed(; - dim = 2, - Ptype = Float64, - tol = 1.0, - C = 0.2, - K = 1.0, - iterations = 100, - initialpos = NL.Point{dim, Ptype}[], - seed = 1, - fixed = false, -) - if !isempty(initialpos) - initialpos = NL.Point.(initialpos) - Ptype = eltype(eltype(initialpos)) - # TODO fix initial pos if list has NL.Points of multiple types - Ptype == Any && error("Please provide list of NL.Point{N,T} with same T") - dim = length(eltype(initialpos)) - end - return SFDP_Fixed{dim, Ptype, typeof(tol)}( - tol, - C, - K, - iterations, - initialpos, - seed, - fixed, - ) -end - -function Base.iterate( - iter::NL.LayoutIterator{SFDP_Fixed{Dim, Ptype, T}}, -) where {Dim, Ptype, T} - algo, adj_matrix = iter.algorithm, iter.adj_matrix - N = size(adj_matrix, 1) - M = length(algo.initialpos) - rng = NL.Random.MersenneTwister(algo.seed) - startpos = Vector{NL.Point{Dim, Ptype}}(undef, N) - startposbounds = - (min = zero(NL.Point{Dim, Ptype}), max = zero((NL.Point{Dim, Ptype})) .+ one(Ptype)) - # take the first - for i in 1:min(N, M) - startpos[i] = algo.initialpos[i] - end - - # create bounds for random initial positions - if M > 0 - startposbounds = ( - min = NL.Point{Dim, Ptype}([ - minimum((p[d] for p in startpos[1:min(N, M)])) for d in 1:Dim - ]), - max = NL.Point{Dim, Ptype}([ - maximum((p[d] for p in startpos[1:min(N, M)])) for d in 1:Dim - ]), - ) - end - - # fill the rest with random NL.Points - for i in (M + 1):N - startpos[i] = - (startposbounds.max .- startposbounds.min) .* - (2 .* rand(rng, NL.Point{Dim, Ptype}) .- 1) .+ startposbounds.min - end - # iteratorstate: (#iter, energy, step, progress, old pos, stopflag) - return startpos, (1, typemax(T), one(T), 0, startpos, false) -end - -function Base.iterate(iter::NL.LayoutIterator{<:SFDP_Fixed}, state) - algo, adj_matrix = iter.algorithm, iter.adj_matrix - iter, energy0, step, progress, locs0, stopflag = state - K, C, tol = algo.K, algo.C, algo.tol - - # stop if stopflag (tol reached) or nr of iterations reached - if iter >= algo.iterations || stopflag - return nothing - end - - locs = copy(locs0) - energy = zero(energy0) - Ftype = eltype(locs) - N = size(adj_matrix, 1) - M = algo.fixed ? length(algo.initialpos) : 0 - for i in 1:N - force = zero(Ftype) - for j in 1:N - i == j && continue - if adj_matrix[i, j] == 1 - # Attractive forces for adjacent nodes - force += Ftype( - f_attr(locs[i], locs[j], K) .* - ((locs[j] .- locs[i]) / norm(locs[j] .- locs[i])), - ) - else - # Repulsive forces - force += Ftype( - f_repln(locs[i], locs[j], C, K) .* - ((locs[j] .- locs[i]) / norm(locs[j] .- locs[i])), - ) - end - end - if i > M - locs[i] = locs[i] .+ step .* (force ./ norm(force)) - end - energy = energy + norm(force)^2 - end - step, progress = update_step(step, energy, energy0, progress) - - # if the tolerance is reached set stopflag to keep claculated NL.Point but stop next iteration - if dist_tolerance(locs, locs0, K, tol) - stopflag = true - end - - return locs, (iter + 1, energy, step, progress, locs, stopflag) -end - -# Calculate Attractive force -f_attr(a, b, K) = (norm(a .- b) .^ 2) ./ K -# Calculate Repulsive force -f_repln(a, b, C, K) = -C .* (K^2) / norm(a .- b) - -function update_step(step, energy::T, energy0, progress) where {T} - # cooldown step - t = T(0.9) - if energy < energy0 - progress = progress + 1 - if progress >= 5 - progress = 0 - step = step / t - end - else - progress = 0 - step = t * step - end - return step, progress -end - -function dist_tolerance(locs, locs0, K, tol) - # check whether the layout is optimal - for i in 1:size(locs, 1) - if norm(locs[i] .- locs0[i]) >= K * tol - return false - end - end - return true -end diff --git a/test/Project.toml b/test/Project.toml index 02e2b7e..4d30cea 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" PowerSystems = "bcd98974-b02a-5e2f-9ee0-a103f5c450dd" PowerSystemsMaps = "e146f72c-d4f8-45d3-b3ca-12a16cb68a38" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/test/runtests.jl b/test/runtests.jl index f4987b7..6a2bf5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,11 +23,13 @@ LOG_LEVELS = Dict( function make_test_sys() sys = System(joinpath(TEST_DIR, "test_data", "case_ACTIVSg200.m")) - locs = CSV.read(joinpath(TEST_DIR, "test_data", "bus_locs.csv"), DataFrames.DataFrame) - for row in DataFrames.eachrow(locs) + locs = PSY.CSV.read( + joinpath(TEST_DIR, "test_data", "bus_locs.csv"), + PSY.DataFrames.DataFrame, + ) + for row in PSY.DataFrames.eachrow(locs) bus = first(get_components(x -> occursin(row.name, get_name(x)), Bus, sys)) if !isnothing(bus) - #set_ext!(bus, Dict("latitude" => row.latitude, "longitude" => row.longitude)) add_supplemental_attribute!( sys, bus, diff --git a/test/test_maps.jl b/test/test_maps.jl index 7a88b5c..f79babf 100644 --- a/test/test_maps.jl +++ b/test/test_maps.jl @@ -4,11 +4,8 @@ sys = make_test_sys() g = make_graph(sys; K = 0.01) @test typeof(g) <: PSM.MetaGraphsNext.MetaGraph @test length(PSM.get_prop(g, :x)) == 200 - shp = PSM.Shapefile.shapes( - PSM.Shapefile.Table( - joinpath(TEST_DIR, "test_data", "IL_BNDY_County", "IL_BNDY_County_Py.shp"), - ), - ) + shapefile = joinpath(TEST_DIR, "test_data", "IL_BNDY_County", "IL_BNDY_County_Py.shp") + shp = PSM.Shapefile.shapes(PSM.Shapefile.Table(shapefile)) shp = PSM.lonlat_to_webmercator(shp) #adjust coordinates @test length(shp) == 102 @@ -39,4 +36,25 @@ sys = make_test_sys() buffer = 0.4e4, ) @test typeof(p) == PSM.Plots.Plot{PSM.Plots.GRBackend} + + p = plot_map( + sys, + joinpath(TEST_DIR, "test_data", "IL_BNDY_County", "IL_BNDY_County_Py.shp"); + map_fillcolor = "grey", + map_background_color = "white", + map_linecolor = "darkgrey", + map_axis = nothing, + map_border = :none, + map_label = "", + map_legend_font_color = :red, + nodesize = 3.0, + linecolor = "blue", + linewidth = 0.6, + lines = true, + nodealpha = 1.0, + shownodelegend = true, + size = (1500, 800), + buffer = 0.4e4, + ) + @test typeof(p) == PSM.Plots.Plot{PSM.Plots.GRBackend} end