MakieUSPop.jl

# setup
# import Pkg;
# Pkg.pkg"add Makie#master PlotUtils GeoInterface GeoJSON"

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")

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...)
    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)