Merge branch 'jpsl/geojson-0.6' of https://github.com/jeremiahpslewis/GeoMakie.jl into jeremiahpslewis-jpsl/geojson-0.6

Author: SimonDanisch

Project.toml

@@ -21,7 +21,7 @@ StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 Colors = "0.12"
 Downloads = "1"
 GeoInterface = "0.5, 1.0"
-GeoJSON = "0.5"
+GeoJSON = "0.6"
 GeometryBasics = "0.4.1"
 ImageIO = "0.6"
 Makie = "0.17.1"

docs/src/index.md

@@ -15,11 +15,10 @@ You can install it from the REPL like so:
 ## GeoAxis
 Using GeoMakie.jl is straightforward, although it does assume basic knowledge of the Makie.jl ecosystem.
 
-GeoMakie.jl provides an axis for plotting geospatial data, [`GeoAxis`](@ref), and also the function [`geo2basic`](@ref) that converts an output of GeoJSON to a polygon appropriate for plotting. Both are showcased in the examples below.
+GeoMakie.jl provides an axis for plotting geospatial data, [`GeoAxis`](@ref). Both are showcased in the examples below.
 
 ```@docs
 GeoAxis
-geo2basic
 ```
 
 ## Gotchas
@@ -111,7 +110,7 @@ using GeoMakie.GeoJSON
 countries_file = download("https://datahub.io/core/geo-countries/r/countries.geojson")
 countries = GeoJSON.read(read(countries_file, String))
 
-n = length(GeoInterface.features(countries))
+n = length(countries)
 hm = poly!(ax, countries; color= 1:n, colormap = :dense,
     strokecolor = :black, strokewidth = 0.5, overdraw = true,
 )

examples/field_and_countries.jl

@@ -4,10 +4,19 @@
 using Makie, CairoMakie, GeoMakie
 import Downloads
 using GeoMakie.GeoJSON
+using GeometryBasics
+using GeoInterface
 
 # https://datahub.io/core/geo-countries#curl # download data from here
 worldCountries = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json"), String))
-n = length(GeoMakie.GeoInterface.features(worldCountries))
+worldCountries = GeoMakie.geoJSONtraitParse(worldCountries)
+
+# TODO: Type issue, Makie won't plot mixed vector of polygons and multipolygons
+worldCountries_poly = [c for c in worldCountries if GeoInterface.trait(c) == PolygonTrait()]
+worldCountries_multipoly = [c for c in worldCountries if GeoInterface.trait(c) == MultiPolygonTrait()]
+
+n_poly = length(worldCountries_poly)
+n_multipoly = length(worldCountries_multipoly)
 lons = -180:180
 lats = -90:90
 field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]
@@ -25,13 +34,21 @@ ax = GeoAxis(
 hm1 = surface!(ax, lons, lats, field; shading = false)
 
 hm2 = poly!(
-    ax, worldCountries;
-    color= 1:n,
+    ax, worldCountries_poly;
+    color= 1:n_poly,
+    colormap = Reverse(:plasma),
+    strokecolor = :black,
+    strokewidth = 0.25
+)
+
+hm2 = poly!(
+    ax, worldCountries_multipoly;
+    color= (1:n_multipoly) .+ n_poly,
     colormap = Reverse(:plasma),
     strokecolor = :black,
     strokewidth = 0.25
 )
 
-cb = Colorbar(fig[1,2]; colorrange = (1, n), colormap = Reverse(:plasma), label = "variable, color code", height = Relative(0.65))
+cb = Colorbar(fig[1,2]; colorrange = (1, n_poly + n_multipoly), colormap = Reverse(:plasma), label = "variable, color code", height = Relative(0.65))
 
 fig

examples/italy.jl

@@ -1,10 +1,13 @@
 using CairoMakie, GeoMakie
 using GeoMakie.GeoJSON
 using Downloads
+using GeometryBasics
+using GeoInterface
+
 # Acquire data
 it_states = Downloads.download("https://github.com/openpolis/geojson-italy/raw/master/geojson/limits_IT_provinces.geojson")
 geo = GeoJSON.read(read(it_states, String))
-basic = GeoMakie.geo2basic(geo)
+basic = GeoMakie.geoJSONtraitParse(geo)
 
 fig = Figure()
 ga = GeoAxis(fig[1, 1]; dest = "+proj=ortho +lon_0=12.5 +lat_0=42", lonlims=(12, 13), latlims = (40, 44))

examples/old/tick_locator_experiment.jl

@@ -66,8 +66,8 @@ function test_tickpad()
         xlimits = (lmin[1], lmax[1])
         ylimits = (lmin[2], lmax[2])
 
-        _xtickvalues, _xticklabels = Makie.MakieLayout.get_ticks(LinearTicks(7), identity, GeoMakie.geoformat_ticklabels, xlimits...)
-        _ytickvalues, _yticklabels = Makie.MakieLayout.get_ticks(LinearTicks(7), identity, GeoMakie.geoformat_ticklabels, ylimits...)
+        _xtickvalues, _xticklabels = Makie.get_ticks(LinearTicks(7), identity, GeoMakie.geoformat_ticklabels, xlimits...)
+        _ytickvalues, _yticklabels = Makie.get_ticks(LinearTicks(7), identity, GeoMakie.geoformat_ticklabels, ylimits...)
 
         _xtickpos_in_inputspace = Point2f.(_xtickvalues, ylimits[1])
         _ytickpos_in_inputspace = Point2f.(xlimits[1], _ytickvalues)

examples/specialized/graph_on_usa.jl

@@ -8,7 +8,7 @@ states_geo = GeoJSON.read(read(states, String))
 # Get rid of Alaska
 filter!((x->!(x.properties["name"] ∈ ("Alaska", "Hawaii", "Puerto Rico"))), states_geo.features)
 
-n = length(GeoInterface.features(states_geo))
+n = length(states_geo)
 
 #g = wheel_graph(10)
 gpos = Dict(

examples/world_population.jl

@@ -1,6 +1,7 @@
 # This example was contributed by Martijn Visser (@visr)
 using Makie, CairoMakie, GeoMakie
 using GeoMakie.GeoJSON
+using GeometryBasics
 using Downloads
 
 source = "+proj=longlat +datum=WGS84"
@@ -19,10 +20,12 @@ ga.yticklabelsvisible[] = false
 url = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/"
 land = Downloads.download(url * "ne_110m_land.geojson")
 land_geo = GeoJSON.read(read(land, String))
+land_geo = GeoMakie.geoJSONtraitParse(land_geo)
 poly!(ga, land_geo, color=:black)
 
 pop = Downloads.download(url * "ne_10m_populated_places_simple.geojson")
 pop_geo = GeoJSON.read(read(pop, String))
-scatter!(ga, GeoMakie.geo2basic(pop_geo), color="lightgrey", markersize=1.2)
+pop_geo = GeoMakie.geoJSONtraitParse(pop_geo)
+scatter!(ga, pop_geo, color="lightgrey", markersize=1.2)
 
 fig

src/GeoMakie.jl

@@ -9,12 +9,13 @@ using GeometryBasics, Colors, ImageIO
 using Makie
 
 import Makie: convert_arguments, convert_attribute, to_value, automatic
-using Makie.MakieLayout, Makie.FileIO, Makie.MakieLayout.GridLayoutBase, Makie.DocStringExtensions
-using Makie.MakieLayout: Formatting
-using Makie.MakieLayout.GridLayoutBase: Side
+using Makie, Makie.FileIO, Makie.GridLayoutBase, Makie.DocStringExtensions
+using Makie: Formatting
+using Makie.GridLayoutBase: Side
 
 using GeoJSON
-using GeoInterface: GeoInterface, coordinates, AbstractPolygon, AbstractMultiPolygon, features, geometry
+using GeoInterface: GeoInterface, coordinates, getfeature
+using GeometryBasics: Polygon, MultiPolygon
 export GeoInterface
 
 
@@ -33,6 +34,7 @@ Makie.to_colormap(::Nothing) = nothing
 # Resolve import conflicts
 import Makie: rotate! # use LinearAlgebra.rotate! otherwise
 
+include("patch.jl")# TODO: Remove this!
 include("conversions.jl")
 include("data.jl")
 include("utils.jl")
@@ -44,6 +46,6 @@ export FileIO
 
 include("geoaxis.jl")
 
-export GeoAxis, geo2basic, datalims, datalims!, automatic
+export GeoAxis, datalims, datalims!, automatic
 
 end # module

src/conversions.jl

@@ -1,62 +1,10 @@
 # # Helper functions
 to_point2(a::Vector{<: Number}) = Point2f(a[1], a[2])
 
-"""
-    geo2basic(obj::GeoInterface)
-
-Converts any GeoInterface object to the equivalent GeometryBasics, which can be plotted by Makie.
-"""
-function geo2basic(poly::GeoInterface.AbstractPolygon)
-    polygon_coordinates = GeoInterface.coordinates(poly)
-    linestrings = map(x-> to_point2.(x), polygon_coordinates)
-    return GeometryBasics.Polygon(linestrings[1], linestrings[2:end])
-end
-
-function geo2basic(poly::Vector{Vector{Vector{Float64}}})
-    linestrings = map(x-> to_point2.(x), poly)
-    return GeometryBasics.Polygon(linestrings[1], linestrings[2:end])
-end
-
-function geo2basic(linestring::Vector{Vector{Float64}})
-    return to_point2.(linestring)
-end
-
-function geo2basic(point::GeoInterface.AbstractPoint)
-    return to_point2(GeoInterface.coordinates(point))
-end
-
-function geo2basic(multi_poly::GeoInterface.AbstractMultiPolygon)
-    polygons = GeoInterface.coordinates(multi_poly)
-    return geo2basic.(polygons)
-end
-
-geo2basic(feature::GeoInterface.Feature) = geo2basic(GeoInterface.geometry(feature))
-function geo2basic(feature::GeoInterface.AbstractLineString)
-    return GeometryBasics.LineString(geo2basic(GeoInterface.coordinates(feature)))
-end
-
-to_multipoly(poly::Polygon) = GeometryBasics.MultiPolygon([poly])
-to_multipoly(mp::MultiPolygon) = mp
-to_multipoly(any) = GeometryBasics.MultiPolygon(any)
-
-# Only converts polygons and multipolygons
-function geo2basic(fc::GeoInterface.AbstractFeatureCollection)
-    return map(geo2basic, GeoInterface.features(fc))
-end
-
-function convert_arguments(P::Type{<: Union{Poly, Mesh}}, geom::GeoInterface.AbstractGeometry)
-    return convert_arguments(P, geo2basic(geom))
-end
-
-function convert_arguments(P::Type{<:Poly}, geom::GeoInterface.AbstractFeatureCollection)
-    return convert_arguments(P, to_multipoly.(geo2basic(geom)))
-end
-
-
 # set the default plot type for Vectors of polygons,
 # so that they are plotted using the most efficient method!
-plottype(::Vector{<: GeoInterface.AbstractMultiPolygon}) = Mesh
-plottype(::Vector{<: GeoInterface.AbstractPolygon}) = Mesh
+plottype(::Vector{<: GeometryBasics.MultiPolygon}) = Mesh
+plottype(::Vector{<: GeometryBasics.Polygon}) = Mesh
 
 # Define a specialized conversion function for images with intervals
 # This means that one can plot images with intervals into GeoMakie

src/data.jl

@@ -11,14 +11,14 @@ end
 function coastlines()
     return get!(LOAD_CACHE, "coastlines") do
         geometry = GeoJSON.read(read(assetpath("vector", "110m_coastline.geojson"), String))
-        return geo2basic(geometry)
+        return geometry
     end
 end
 
 
 function land()
     return get!(LOAD_CACHE, "land") do
         geometry = GeoJSON.read(read(assetpath("vector", "110m_land.geojson"), String))
-        return geo2basic(geometry)
+        return geometry
     end
 end

src/geoaxis.jl

@@ -1,5 +1,5 @@
 using Makie: left, right, top, bottom
-using Makie.MakieLayout: height, width
+using Makie: height, width
 
 """
     GeoAxis(fig_or_scene; kwargs...) → ax::Axis
@@ -121,8 +121,7 @@ function GeoAxis(args...;
     Makie.Observables.connect!(ax.scene.transformation.transform_func, _transformation)
 
     # Plot coastlines
-    coast_line = Makie.convert_arguments(PointBased(), GeoMakie.coastlines())[1]
-
+    coast_line = GeoMakie.geoJSONtraitParse(GeoMakie.coastlines())
     coastplot = lines!(ax, coast_line; color = :black, coastline_attributes...)
     translate!(coastplot, 0, 0, 99) # ensure they are on top of other plotted elements
     xprot = ax.xaxis.protrusion[]
@@ -204,11 +203,11 @@ function draw_geoticks!(ax::Axis, hijacked_observables, line_density, remove_ove
         xlimits[] = (lmin[1], lmax[1])
         ylimits[] = (lmin[2], lmax[2])
 
-        _xtickvalues, _xticklabels = Makie.MakieLayout.get_ticks(xticks, identity, xtickformat, xlimits[]...)
-        _ytickvalues, _yticklabels = Makie.MakieLayout.get_ticks(yticks, identity, ytickformat, ylimits[]...)
+        _xtickvalues, _xticklabels = Makie.get_ticks(xticks, identity, xtickformat, xlimits[]...)
+        _ytickvalues, _yticklabels = Makie.get_ticks(yticks, identity, ytickformat, ylimits[]...)
 
-        _xminortickvalues = Makie.MakieLayout.get_minor_tickvalues(xminor, identity, _xtickvalues, xlimits[]...)
-        _yminortickvalues = Makie.MakieLayout.get_minor_tickvalues(yminor, identity, _ytickvalues, ylimits[]...)
+        _xminortickvalues = Makie.get_minor_tickvalues(xminor, identity, _xtickvalues, xlimits[]...)
+        _yminortickvalues = Makie.get_minor_tickvalues(yminor, identity, _ytickvalues, ylimits[]...)
 
         _xtickpos_in_inputspace = Point2f.(_xtickvalues, ylimits[][1])
         _ytickpos_in_inputspace = Point2f.(xlimits[][1], _ytickvalues)

src/patch.jl

@@ -0,0 +1,40 @@
+using GeoInterface
+using GeometryBasics
+
+# Coerce all coordinates to Floats
+#### Hack for type inference issues for integers, related to ####################################
+#### FIX FOR https://github.com/JuliaGeometry/GeometryBasics.jl/issues/142 #####################
+#### Should be fixed by https://github.com/JuliaGeometry/GeometryBasics.jl/pull/173, perhaps. ##
+function GeoInterface.convert(::Type{LineString}, type::LineStringTrait, geom)
+    dim = Int(ncoord(geom))
+    return LineString([Point{dim}(GeoInterface.coordinates(p) * 1.0) for p in getgeom(geom)])
+end
+
+function GeoInterface.convert(::Type{Point}, type::PointTrait, geom)
+    dim = Int(ncoord(geom))
+    return Point{dim}(GeoInterface.coordinates(geom) * 1.0)
+end
+### Hack over...
+
+trait_type_pairs = [PolygonTrait() Polygon;
+                    MultiPolygonTrait() MultiPolygon;
+                    LineStringTrait() LineString;
+                    PointTrait() Point]
+
+# Convert all geoJSON objects to the appropriate GeometryBasics type based on trait
+function geoJSONtraitParse(feature::GeoJSON.Feature)
+    geometry = GeoInterface.geometry(feature)
+    trait_matches = trait_type_pairs[:,1] .== (GeoInterface.trait(geometry),)
+
+    if sum(trait_matches) != 1
+        @warn "GeoMakie.geoJSONtraitParse: Unknown geometry type $(GeoInterface.trait(geometry))"
+        return nothing
+    end
+
+    return GeoInterface.convert(trait_type_pairs[trait_matches, 2][1], geometry)
+end
+
+function geoJSONtraitParse(featureCollection::GeoJSON.FeatureCollection)
+    Vector{}
+    return [geoJSONtraitParse(f) for f in featureCollection]
+end

src/utils.jl

@@ -198,7 +198,7 @@ function latitude_format(nums)
 end
 
 function _replace_if_automatic(typ::Type{T}, attribute::Symbol, auto) where T
-    default_attr_vals = Makie.MakieLayout.default_attribute_values(T, nothing)
+    default_attr_vals = Makie.default_attribute_values(T, nothing)
 
     if to_value(get(default_attr_vals, attribute, automatic)) == automatic
         return auto

test/runtests.jl

@@ -1,4 +1,4 @@
-using GeoMakie, CairoMakie, Test
+using GeoMakie, GeometryBasics, CairoMakie, Test
 
 Makie.set_theme!(Theme(
     Heatmap = (rasterize = 5,),
@@ -18,6 +18,11 @@ Makie.set_theme!(Theme(
     # display(fig)
 end
 
+@testset "geoJSONtraitParse" begin
+    @test GeoMakie.geoJSONtraitParse(GeoMakie.coastlines()) isa Vector
+    @test GeoMakie.geoJSONtraitParse(GeoMakie.coastlines()[1]) isa GeometryBasics.LineString
+end
+
 @testset "Examples" begin
     geomakie_path = dirname(dirname(pathof(GeoMakie)))
     examples = readdir(joinpath(geomakie_path, "examples"); join = true)