So the Easting values change a bit and the Northing values change drastically. I have tried to recreate the old behaviour but without success so far. What has changed? Is there a way to get the old projection values with the newer versions or do I have to stick to 6.0 and before?
I have a grid file that I want to inversely project with the same settings (it is a rectilinear grid used to model the data after projection). I.e. something like
However this produces an empty grid. The issue appears to be that the region set in -R does not contain the coordinates used in the projection. This is because the dataset is quite small in scale (a few kilometers) and when I initially selected the projection I simply rounded to the nearest lat/lon even though they are outside the region of interest. I thought with these small distances it should not be an issue. Is this problem a bug or by design? Is there a way to fix this other than extending the region, creating a large grid (which will be enormous because the increment is 100m) and then cutting down again.
When setting the region in grdproject is there a convenient/semi-automatic way to set the grid parameters (-R -D) so they are as close as possible to the extents of my grid after projection. I know there is -Rgrid, but as I understand it that would have to have the information after projection (i.e. boundaries and increment in lat/lon). So at the moment I have to determine the projections of the corners manually and fill those in the command line/script. For my purposes it would be ideal if I could specify the grid, grdproject would determine the extents and use the same number of divisions as in my grid. This should hopefully result in a projected grid that represents the input well.
Hope all this makes any sense, happy to answer any questions that arise.
I do recall we addressed a bug in the Mercator projection for when a specific lon/lat center was given. Since you are setting (-81/46) as the center of the projection, you are now getting “northings” relative to 46, not the Equator. Most folks projecting data to easting and northing do UTM, not Mercator, so we had not really implemented that very carefully for regular Mercator. I will have to do some more checking so I understand exactly what we changed to answer your question. If you use -Jm-18/0/1:1 you will see your northings are closer to what you had before.but still quite different.
However, using Mercator and terms like “northing” is not good since these are not distances in the normal sense of those terms. Recommend you look into UTM.
Thanks for those quick replies, using the proj4 syntax does give the previous results with gmt 6.1 and mapproject.
Regarding grdproject, using it with proj4 syntax now gives the expected while using -Jm-81/0/1:1 still gives an empty. This is despite the fact that grdproject with proj4 syntax claims to convert it to -Jm-81/0/1:1 and appears to understand the coordinates in the same way. I put both outputs below, so you can maybe see what is going on.
/usr/local/gmt-6.1.0/build/src/gmt grdproject -Jm-81/0/1:1 -n+c -D100+ -R-81.6/-81.4/46.4/46.6 -F -C -I tmp.nc?Conductivity[0] -Gmerge2.grd -V
grdproject [INFORMATION]: Processing input grid
grdproject [INFORMATION]: No range attribute, guessing registration to be pixel
grdproject [INFORMATION]: netCDF grid tmp.nc information has zmin = zmax = NaN. Reset to 0/0.
grdproject [INFORMATION]: Cartesian input grid
grdproject [INFORMATION]: Reading grid from file tmp.nc?Conductivity[0]
grdproject [INFORMATION]: Cartesian input grid
grdproject [INFORMATION]: Set boundary condition for all edges: natural
grdproject [INFORMATION]: Set boundary condition for left edge: natural
grdproject [INFORMATION]: Set boundary condition for right edge: natural
grdproject [INFORMATION]: Set boundary condition for bottom edge: natural
grdproject [INFORMATION]: Set boundary condition for top edge: natural
grdproject [INFORMATION]: Given n_columns implies x_inc = 0.002
grdproject [INFORMATION]: Given n_rows implies y_inc = 0.002
grdproject [INFORMATION]: Grid projection from size 0x0 to 100x100
grdproject [INFORMATION]: Transform (-81.6/-81.4/46.4/46.6) <-- (-66791.694476/-44527.7963173/5813726.13738/5845966.88134) [m]
grdproject [INFORMATION]: Writing grid to file merge2.grd
grdproject [INFORMATION]: Proj4 string to be converted to WKT:
+proj=longlat +ellps=WGS84 +no_defs
grdproject [INFORMATION]: WKT converted from proj4 string:
GEOGCS[“WGS 84”,
DATUM[“unknown”,
SPHEROID[“WGS84”,6378137,298.257223563]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]]
grdproject [INFORMATION]: No valid values in grid [merge2.grd]
/usr/local/gmt-6.1.0/build/src/gmt grdproject -J+proj=merc+lat_ts=46+lon_0=-81 -n+c -D100+ -R-81.6/-81.4/46.4/46.6 -F -C -I tmp.nc?Conductivity[0] -Gmerge.grd -V
grdproject [INFORMATION]: Converted to -J syntax = -Jm-81/0/1:1
grdproject [INFORMATION]: Processing input grid
grdproject [INFORMATION]: No range attribute, guessing registration to be pixel
grdproject [INFORMATION]: netCDF grid tmp.nc information has zmin = zmax = NaN. Reset to 0/0.
grdproject [INFORMATION]: Cartesian input grid
grdproject [INFORMATION]: Reading grid from file tmp.nc?Conductivity[0]
grdproject [INFORMATION]: Cartesian input grid
grdproject [INFORMATION]: Set boundary condition for all edges: natural
grdproject [INFORMATION]: Set boundary condition for left edge: natural
grdproject [INFORMATION]: Set boundary condition for right edge: natural
grdproject [INFORMATION]: Set boundary condition for bottom edge: natural
grdproject [INFORMATION]: Set boundary condition for top edge: natural
grdproject [INFORMATION]: Given n_columns implies x_inc = 0.002
grdproject [INFORMATION]: Given n_rows implies y_inc = 0.002
grdproject [INFORMATION]: Grid projection from size 0x0 to 100x100
grdproject [INFORMATION]: Transform (-81.6/-81.4/46.4/46.6) <-- (-66791.694476/-44527.7963173/5813726.13738/5845966.88134) [m]
grdproject [INFORMATION]: Writing grid to file merge.grd
grdproject [INFORMATION]: Proj4 string to be converted to WKT:
+proj=longlat +ellps=WGS84 +no_defs
grdproject [INFORMATION]: WKT converted from proj4 string:
GEOGCS[“WGS 84”,
DATUM[“unknown”,
SPHEROID[“WGS84”,6378137,298.257223563]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433]]
Might be a little dubious but it doesn’t claim to use the -J syntax, only that it was able to convert the proj4 string to it. In fact when using -J+proj4 (or EPSG) the conversions are done with the proj4 library, not the GMT projection machinery.
So (and I’m hearing you thinking out loud) why not use proj4 all the time? Well, because there are lots of ramifications that need to be carried on in that transition and they are not all done (creating the frames in maps is probably the harder one).