Plotting countries with GeoMakie: reduce size and more

filchristou · April 26, 2023, 7:52pm

I am trying visualize some specific countries with GeoMakie.jl.
I 've seen Examples · GeoMakie.jl and I tried to do something similar:

begin
	# Acquire data
	ger_states = Downloads.download("https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/4_niedrig.geo.json")
	geoger = GeoJSON.read(read(ger_states, String))

	fra_states = Downloads.download("https://raw.githubusercontent.com/gregoiredavid/france-geojson/master/regions.geojson")
	geofra = GeoJSON.read(read(fra_states, String))
	
	fig = Figure()
	ga = GeoAxis(fig[1, 1]; source = "+proj=longlat +datum=WGS84", dest = "+proj=lcc +lon_0=2 +lat_1=10 lat_2=20", lonlims=(12, 13), latlims = (30, 44), title="Simulation topology", coastlines=true)
	poly!(ga, geoger; strokewidth = 1, color=:lightblue)
	poly!(ga, geofra; strokewidth = 1, color=:royalblue)
	datalims!(ga)

	save("gerfratop.pdf", fig)
	fig
end

It would be great if you could help me refine this figure a bit. Following are my problems:

My biggest problem is that this figure is 43MB. As a result, it’s not very usable. Is there a way to reduce the size to the KB range ? Maybe somehow sample the found geometry in the json files ?
The details of the coastlines are too coarse. How can I increase it to make it much the shapefiles used ? Also the coastlines shouldn’t escape the frame.
How can I paint the sea ?
Bottom right the latitude text is cut in half and it’s a bit ugly.
How can I also add lakes, and a sketch of the neighboring countries borders (without the inner states). Do I need to find the corresponding shapefile for each one of them, or can it be automated (similar to the coastilines, which uses the Natural Earth dataset) ? I tried doing

worldCountries = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/johan/world.geo.json/master/countries.geo.json"), String))
poly!(ga, worldCountries; color= :grey)

…but it didn’t ended up as desired:

The style I am aiming to achieve is something similar to the following but for both France and Germany (ignore the overlaid graph).

Thanks a lot in advance !

Reproducibility

Julia v1.8.5 and GeoMakie 0.5.0

joa-quim · April 26, 2023, 8:24pm

If you want to try GMT, you can:
Note image is a bit ugly due to another ghostscript bug in version 10 (versions 9.xx work fine). Using pdf as output workaround the issue. png file size ~500 KB.
(EDIT: high quality pdf = 263 KB)

using GMT

julia> ger_states = gmtread("https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/4_niedrig.geo.json");
        This file has islands (holes in polygons).
        Use `gmtread(..., no_islands=true)` to ignore them.

julia> fra_states = gmtread("https://raw.githubusercontent.com/gregoiredavid/france-geojson/master/regions.geojson");
        This file has islands (holes in polygons).
        Use `gmtread(..., no_islands=true)` to ignore them.

julia> coast(proj="+proj=lcc +lon_0=2 +lat_1=10 +lat_2=20", region=(-5,15,41,55), coast=0.5, water=:lightgreen)
julia> plot!(ger_states, lc=:black, fill=:blue)
julia> plot!(fra_states, lc=:black, fill=:lightblue, show=1)

asinghvi17 · April 26, 2023, 10:53pm

Yeah, GMT’s land-sea mask probably has a lot of the details you’re aiming for.

From a GeoMakie perspective, it’s going through a major refactor right now, so things should change in a couple of weeks or so.

To get the file size to be smaller, set the rasterize attribute in your poly plots. rasterize = 5 or so should be pretty smooth, and will basically render the polygons to an image and then plot that image instead of plotting the polygons directly. This is theoretically not “true vector” but is good for most usecases, since printers can’t get that level of accuracy anyway.
To get smoother coastlines, you can use natural-earth-vector/ne_10m_coastline.geojson at master · nvkelso/natural-earth-vector · GitHub from Natural Earth. You would have to plot this yourself though. I am working on getting more natural earth data into GeoMakie but that will take some time, see e.g. https://github.com/NaturalEarth.jl
You could use the land-sea-river mask from GitHub - JuliaGeo/GeoDatasets.jl: Common geographic data sets to paint the sea.
In general the tick label spacing was always a bit buggy. We came up with a nice method to figure out tick label spacing at MakieCon, so I will try and get a patch out which uses that soonish.
For lakes, you can use natural-earth-vector/ne_10m_lakes_europe.geojson at master · nvkelso/natural-earth-vector · GitHub. For neighbouring countries, your current approach works. However, when plotting the country outlines, set xautolimits=false, yautolimits=false. This will make sure that those plots aren’t taken into account when checking the limits!

Hope these help. I will try to get a patch out for the tick labels soonish.

joa-quim · April 26, 2023, 11:03pm

Just to clarify this, it’s not a mask. It’s a pure vector operation (polygon painting in postscript). That is why the PDF version is so small ~260 KB)

jules · April 27, 2023, 5:52am

Are the polygons rendering as meshes again in this example? Otherwise Makie should also just export them as plain polys in pdf files, keeping the size small. Some types revert to mesh drawing though, maybe this is happening here. That blows up the filesize and can be fixed in cairomakie.

asinghvi17 · April 28, 2023, 1:06pm

Definitely possible…@filchristou, which version of CairoMakie are you using? You can find out by:

using Pkg
Pkg.status("CairoMakie")

filchristou · April 28, 2023, 3:37pm

Thanks for the answers. I will work my way through them and I will update you!

I am using CairoMakie 0.10.4

filchristou · April 28, 2023, 6:50pm

@asinghvi17

rasterize=5 appears to work okey; Now I am not sure how to check what @jules is saying.
I tried it but I am getting some errors with axisautolimits=false similar to this issue I made {x,y}autolimits=false breaks example · Issue #166 · MakieOrg/GeoMakie.jl · GitHub . Also plotting is quite slow. Can I somehow filter only the region I am interested in, so that I don’t wait for the plotting of the whole earth ?
The following appears to work, but again quite slow for high quality.
Also I would prefer a solution with poly! as it’s messy handling it as a field.
ok. I will completely remove this information. However, I couldn’t remove the grid and I made an issue make grid invisible · Issue #167 · MakieOrg/GeoMakie.jl · GitHub
Lakes appears to be located wrong. I made an issue wrong coordinate transformation ? German lakes · Issue #168 · MakieOrg/GeoMakie.jl · GitHub

So what I have now is the following:

begin
	allcountries = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_10m_admin_0_countries.geojson"), String))

	geoger = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/4_niedrig.geo.json"), String))
	
	geofra = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/gregoiredavid/france-geojson/master/regions.geojson"), String))

	lakes = GeoJSON.read(read(Downloads.download("https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/ne_10m_lakes_europe.geojson"), String))

	lontemp,lantemp,dattmep = GeoDatasets.landseamask(;resolution='l',grid=5)
	clamp!(dattmep, 0x00, 0x01)
end

mymapfig = let
	fig = Figure()
	
	ga = GeoAxis(fig[1, 1]; dest = "+proj=merc", lonlims=(12, 13), latlims = (30, 44), title="Simulation topology", xgridvisible=false, ygridvisible=false, xticklabelsvisible=false, yticklabelsvisible=false, xticksvisible=false, yticksvisible=false)

	contourf!(ga, lontemp, lantemp, dattmep; colormap=ColorScheme(range(colorant"lightblue", colorant"lightgrey", length=2)),  xautolimits=false, yautolimits=false)
	poly!(ga, allcountries; color= :grey, xautolimits=false, yautolimits=false, rasterize=5)
	poly!(ga, geoger; strokewidth = 1, color=:gold, rasterize=5)
	poly!(ga, geofra; strokewidth = 1, color=:peru, rasterize=5)
	poly!(ga, lakes; strokewidth = 1, color=:blue, rasterize = 5,  xautolimits=false, yautolimits=false)
	
	datalims!(ga)
	
	fig
end

it produces this 2.7MB fig:

Remarks:
i) lakes are positioned wrongly
ii) grind lines are visible
iii) there is still some white between the gray countries and the blue see. I guess this is due to using different datasets.
iv) I still need land boundaries. I didn’t find a geojson file in the natural earth repo. However there are many boundaries here with formats like *.cpg, *.dbf, *.prj, *.shp, *.shx. Is there a way to read any of them and use them with GeoMakie ?
v) generation is slow. I think there are some low-hanging fruits by filtering the coordinates only to the frame we are interested into (?)

@joa-quim thanks for the solution. GMT looks powerful and mature, but I need GeoMakie here because on top I intend to use some GraphMakie recipes.

joa-quim · April 28, 2023, 7:01pm

Sure, it’s your choice. But I can’t understand how the above even works because this

is completely wrong (and I don’t think the figure is in Mercator projection).

filchristou · April 29, 2023, 7:55am

I think because I {x,y}autolimits are automatically set to true, GeoMakie adapts the axis for both german and french topology. However I need to supploy something like lonlims=(12, 13), latlims = (30, 44) in the signature because it’s a required argument.

asinghvi17 · May 17, 2023, 8:14pm

Hmm…I got this kinda working, it looks like the huge file size is a result of a wrong dispatch in CairoMakie

The generated PDF is 430 KB, so not bad. This is on an experimental branch and has quite a few other issues as well, but I just wanted to check that the aspect ratio seems correct? If so I can actually try and update the master version of geoaxis with this.