Closed consumere closed 2 months ago
I can replicate that. For UTM projections you probably want to restrict the longitudes to that particular zone though? And you can restrict latitudes to (-90, 90).
The long term solution to this is to define bounding boxes for each projection, individually.
Yes, i was wondering, if it wouldnt be much effort to implemnt this kind of transformations, since GeoMakie already uses proj strings.
You are already using those, no? It's just a matter of passing the correct limits to GeoAxis then.
No, i would expect a grid with X and Y coordinates. Like this:
import ArchGDAL
srccrs = "+proj=longlat +datum=WGS84 +no_defs"
destcrs = "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
xc,yc = 11.5, 50.5
pt = ArchGDAL.createpoint(xc,yc)
#Geometry: POINT (11.5 50.5)
ArchGDAL.reproject(pt,
ProjString(srccrs),
ProjString(destcrs)
)
#Geometry: POINT (677291.966065446 5597210.83026226)
Ah, for that you can plot in a regular Makie axis! You would have to explicitly reproject your geometries, though.
Got it, thanks... so no internal Raster transformation. for example. this works.
r isa Raster
#true
#already utm
import Makie
fig = GeoMakie.Figure();
ax = Makie.Axis(fig[1,1],title="UTM 32";)
rsf = GeoMakie.heatmap!(ax,r;)
fig
Yep this works! The idea of GeoMakie is to be as transparent as possible with the least amount of added features - everything that works here except for (a) multiple CRS in a plot and (b) lat-long grids in any projection will work in base Makie. Some code has actually been shifted from here to Makie.jl for that reason.
I tried to transform to a UTM projection.
prj = "+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
fig = GeoMakie.Figure();
ax = GeoMakie.GeoAxis(fig[1,1], dest=prj, title="UTM 32";)
lp = lines!(ax, GeoMakie.coastlines(); transformation = (; translation = (0, 0, 1)))
This is what i get: