I have (multiple) data on a regular latlon grid which I want to project with a Lambert conformal conic projection and interpolate to a regular grid. Is using the Proj4 and Interpolations packages an appropriate option or is there a more tailored alternative?
Would it be possible to provide some more details, some reference or some images illustrating your problem?
I am not an expert in Geodesy, just an interested user, and in the mock-up example below produced with Proj4.jl, I have tried to visualize your problem, probably in a wrong way, but hope it helps getting more feedback.
The left image is a regular Lat-Lon grid in WGS84 and the right image its transformation to North America Lambert conformal conic (LCC) X-Y coordinates.
It does not make immediate sense to define a regular grid on this output. What would be the application?
PS: from @joa-quim’s GMT grdproject reference provided, one may define a regular grid that fits inside the projected LCC irregular grid, then do an inverse projection and perform bicubic spline interpolation (ex: with Dierckx.jl Spline2D) using the available spatial data (ex: temperature T(X,Y)) in the original grid.
Yes, your example illustrates the first part well. The next step is to interpolate the data on the Lambert grid to a grid expanded by the vectors, say,
x = 1e6:1e4:2e6
y = 3e5:1e4:2e6
The aim is to match two datasets with different projections.
If you have a grid in geographical coordinates you can reproject that grid into Lcc and get regular cartesian grid in Lcc meters. Either GMT’s gdproject or GDAL gdalwarp let you do that.
2 Likes
Thanks. I will need to look further into these. Are there any Julia examples of using GMT.grdproject?
For GMT you would need to create a GMTgrid type. I’m sorry but cannot give a pointer to the docs because something screw in the documentation builder and the GMT docs are giving a 404 right now. Hope it recovers soon.
But the GMT’s mat2grid function offers lots of possibilities to create them. It has an online help.
The, assuming G holds the geog grid, the command would be something like (need to complete the parameters for projection using the PROJ4 syntax)
Glamb = grdproject(G, proj=“+proj=lcc+lat_1=…+lat_2=…+lat_0=…”) ;
and eventually even easier if you already have the geog grid as a file on disk. Then the comand would be (assuming that GMT will guess the file format).
Glamg = grdproject(“filename”, proj=“+proj=lcc+lat_1=…+lat_2=…+lat_0=…”) ;
1 Like
@johnbb, if the regular grid is larger than the projected input data, extrapolation will be also required. There will be some geometrical/scale distortion. See also the postscript in my updated post above for additional comments.
The inverse approach is appealing as I have many datasets and the latlon grid is less dense, but covers completely the target area/grid. Proj4 and Interpolations may then be a good alternative.
OK, a full GMT example
# Extract a sub-region from a global grid. Using a low resolution to keep example light
Ggeog = grdcut("@earth_relief_10m", region=(-18,15,32,56));
convert to Lambert
Glamb = grdproject(Ggeog,proj="+proj=lcc+lat_1=39+lat_2=50+lat_0=45+lon_0=0");
Show them. Rafael, you are right about the extrapolation but the trick is to not do it and instead fill those nodes with NaNs.
imshow(Ggeog, shade=true, cmap=:geo, name="geogs.png")
imshow(Glamb, shade=true, cmap=:geo, name="lcc.png")
4 Likes
@joa-quim brilliant plot. Love it.
1 Like