I am trying GeoMakie for a simple visualization of KNMI radar data (Netherlands), which comes in netCDF format. The netCDF provides grid values, the x and y axes and the proj4. I am able to successfully plot that georeferenced with:
ga = GeoAxis(fig[1,1], dest=radar[1].prj, source=radar[1].prj)
heatmap!(ga,radar[1].x, radar[1].y, radar[1].val)
where radar[1].prj is "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378140 +b=6356750 +x_0=0 +y_0=0 +units=km"
It plots a reference grid of latitude and longitude (despite not having a frame around it, which I would prefer).
However, I am struggling to figure out the right way to plot a shapefile of country borders or the built in GeoMakie.coastlines correctly on top of it. Is there a good example available for this?
julia> ga = GeoAxis(fig[1,1], dest=radar[1].prj, source=radar[1].prj)
GeoAxis()
julia> heatmap!(ga,radar[1].x,radar[1].y,radar[1].val)
Heatmap{Tuple{Vector{Float64}, Vector{Float64}, Matrix{Float32}}}
# correctly displays geo-referenced data
julia> countries = Shapefile.Handle(raw"C:\Users\oscar\Documents\Work\Julia\LMA\map\world-administrative-boundaries.shp").shapes
256-element Vector{Union{Missing, Shapefile.Polygon}}:
Polygon(24 Points) # etc...
julia> poly!(ga,countries)
# or:
julia> lines!(ga,GeoMakie.coastlines())
# displays them small, all the way at the north pole
I have tried to play with the source and destination projection strings within the GeoAxis, before/after plotting.
It looks like I cannot set ga.source[] = "+proj=longlat" after plotting the data and before adding the countries, because it also changes how the radar data are interpreted.
Perhaps I should first convert the data grid to latitude longitude (how…). But when I do that, I also wouldn’t need GeoMakie anymore, as projections are not so critical for displaying a limited area up to a few hundred kilometers wide.
I have also tried to load the netCDF with Rasters.jl, as it may have functions for conversion, but it did not pick up the information properly:
Thanks, it is easily converted. This works now. Perhaps it could be worked into an example for the GeoMakie site.
trans = Proj.Transformation(prj,"WGS84",always_xy=true)
xy = [[x[i],y[j]] for i in 1:dims[1], j in 1:dims[2]]
lonlats = trans.(xy)
lon = reshape(getindex.(lonlats,1),dims)
lat = reshape(getindex.(lonlats,2),dims)
fig = nothing; fig = Figure(size = (900,900))
ga = GeoAxis(fig[1,1], limits=((0,10),(49,55)), dest="+proj=merc +lat_ts=52.0")
surface!(ga, lon,lat, values, colorrange=(0,percentile(skipmissing(values),99)))
poly!(ga, countries; strokecolor=(:gray,0.5), strokewidth = 1.25, overdraw = true)
The opposite transformation, from countries or GeoMakie.coastlines() WGS84 to the radar projection seems to need some additional conversion step before:
julia> trans.(GeoMakie.coastlines())
ERROR: ArgumentError: Argument is not a Point geometry
Could anyone help with that?
Update: corrected the alignment of x and y for the coordinate transform. Update 2: Succeeded plotting, here is the result!
