Skip to content

Instantly share code, notes, and snippets.

@benelsen

benelsen/MakieUSPop.jl

Forked from asinghvi17/MakieUSPop.jl
Last active August 3, 2019 11:01
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save benelsen/517ada29777f47564983b9ac047b82d3 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save benelsen/517ada29777f47564983b9ac047b82d3 to your computer and use it in GitHub Desktop.
Download ZIP

Revisions

  1. benelsen revised this gist Aug 3, 2019. 1 changed file with 29 additions and 2 deletions.
    31 changes: 29 additions & 2 deletions MakieUSPop.jl
    View file
    Original file line number Diff line number Diff line change
    @@ -5,13 +5,41 @@
    using Makie
    using GeoInterface, GeoJSON
    using PlotUtils
    using Proj4

    states = download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")

    states_geo = GeoJSON.parse(read(states, String))

    features = states_geo.features

    wgs84 = Projection("+proj=longlat +datum=WGS84 +no_defs")
    albers_us = Projection("+proj=aea +lat_1=29.5 +lat_2=45.5 +lon_0=-96.0 +datum=WGS84 +no_defs")

    function project(geometry::Polygon, from, to)
    newcoords = map(coordinates(geometry)) do ring
    map(ring) do point
    transform(from, to, point)
    end
    end
    Polygon(newcoords)
    end

    function project(geometry::MultiPolygon, from, to)
    newcoords = map(coordinates(geometry)) do element
    map(element) do ring
    map(ring) do point
    transform(from, to, point)
    end
    end
    end
    MultiPolygon(newcoords)
    end

    for feature in features
    feature.geometry = project(geometry(feature), wgs84, albers_us)
    end

    props = getproperty.(features, :properties)

    densities = getindex.(props, "density")
    @@ -27,7 +55,6 @@ function plo(sc::Scene, poly::Polygon; kwargs...)
    end

    function plo(sc, poly::MultiPolygon; kwargs...)
    # map(x -> poly!(sc, Point2f0.(x[1]); kwargs...), poly.coordinates)
    data = map(x -> Point2f0.(x[1]), poly.coordinates)
    poly!(sc, data; kwargs...)
    sc
    @@ -36,7 +63,7 @@ end
    sc = Scene()

    normalize(val::Real; min = cr[1], max = cr[2]) = (val - min) / (max - min)

    for (i, feature) in enumerate(features)
    println(i)
    try
  2. @asinghvi17 asinghvi17 created this gist Aug 3, 2019.
    52 changes: 52 additions & 0 deletions MakieUSPop.jl
    View file
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,52 @@
    # setup
    # import Pkg;
    # Pkg.pkg"add Makie#master PlotUtils GeoInterface GeoJSON"

    using Makie
    using GeoInterface, GeoJSON
    using PlotUtils

    states = download("https://raw.githubusercontent.com/PublicaMundi/MappingAPI/master/data/geojson/us-states.json")

    states_geo = GeoJSON.parse(read(states, String))

    features = states_geo.features

    props = getproperty.(features, :properties)

    densities = getindex.(props, "density")

    cr = let ds = sort(densities); (ds[1], ds[end-1]); end

    Base.convert(poly::GeoInterface.Polygon, ::Type{Point2f0}) = Point2f0.(poly.coordinates[1])

    viridis = cgrad(:viridis; scale = :log10)

    function plo(sc::Scene, poly::Polygon; kwargs...)
    poly!(sc, convert(poly, Point2f0); kwargs...)
    end

    function plo(sc, poly::MultiPolygon; kwargs...)
    # map(x -> poly!(sc, Point2f0.(x[1]); kwargs...), poly.coordinates)
    data = map(x -> Point2f0.(x[1]), poly.coordinates)
    poly!(sc, data; kwargs...)
    sc
    end

    sc = Scene()

    normalize(val::Real; min = cr[1], max = cr[2]) = (val - min) / (max - min)

    for (i, feature) in enumerate(features)
    println(i)
    try
    plo(sc, feature.geometry, color = viridis[normalize(feature.properties["density"])])
    catch err
    @warn("Feature had dimension zero!")
    print(err)
    pop!(sc.plots)
    end
    # display(sc)
    end

    save("usa-pop-density-log.png", sc)