# Mapproject/grdproject issues and changes between versions

**URL:** <https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823>\
**Category:** Q&A\
**Created:** [August 28, 2020, 2:59pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823 "2020-08-28T14:59:12Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![mmoorkamp](https://avatars.discourse-cdn.com/v4/letter/m/f4b2a3/32.png) [@mmoorkamp](https://forum.generic-mapping-tools.org/u/mmoorkamp)\
**Post date:** [August 28, 2020, 2:59pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/1 "2020-08-28T14:59:12Z")

</div>

Hi,

I have been running into a few issues with mapproject/grdproject lately and hope someone here can clarify help.

1. It seems that some defaults have changed with the projections between 6.0.0 and older versions and 6.1.0.

For example this file

-81.50638333 46.47973333 L00S\_0250E\_R\_GT  
-81.51105 46.47666667 L00S\_0250W\_R\_GT  
-81.50403333 46.48143333 L00S\_0500E\_R\_GT  
-81.51328333 46.47516667 L00S\_0500W\_GT  
-81.51573333 46.47348333 L00S\_0750W\_GT  
-81.49945 46.4846 L00S\_1000E\_R\_GT  
-81.518 46.4717 L00S\_1000W\_GT  
-81.49736667 46.48613333 L00S\_1250E\_GT  
-81.52028333 46.47028333 L00S\_1250W\_GT  
-81.52243333 46.46865 L00S\_1500W\_GT

with 6.0 (and older versions as well)

/usr/local/gmt-6.0.0/build/src/gmt mapproject -Jm-81/46/1:1 -F -R-82/-80/45/47 -C tmp | head  
|-39226.1233521|4054500.74761|L00S\_0250E\_R\_GT|  
|-39587.6190061|4054156.88134|L00S\_0250W\_R\_GT|  
|-39044.0845992|4054691.37795|L00S\_0500E\_R\_GT|  
|-39760.6201159|4053988.69261|L00S\_0500W\_GT|  
|-39950.4051987|4053799.9523|L00S\_0750W\_GT|  
|-38689.0447365|4055046.49003|L00S\_1000E\_R\_GT|  
|-40125.9889348|4053600.00726|L00S\_1000W\_GT|  
|-38527.6631216|4055218.44597|L00S\_1250E\_GT|  
|-40302.8632096|4053441.17645|L00S\_1250W\_GT|  
|-40469.4093026|4053258.05983|L00S\_1500W\_GT|

whereas with 6.1 I get

/usr/local/gmt-6.1.0/build/src/gmt mapproject -Jm-81/46/1:1 -F -R-82/-80/45/47 -C tmp | head  
|-39361.7316646|53742.7941827|L00S\_0250E\_R\_GT|  
|-39724.4770423|53397.7391316|L00S\_0250W\_R\_GT|  
|-39179.063587|53934.0835481|L00S\_0500E\_R\_GT|  
|-39898.0762328|53228.9689647|L00S\_0500W\_GT|  
|-40088.5174201|53039.5761576|L00S\_0750W\_GT|  
|-38822.7963189|54290.4232895|L00S\_1000E\_R\_GT|  
|-40264.7081653|52838.9398934|L00S\_1000W\_GT|  
|-38660.8567929|54462.9736875|L00S\_1250E\_GT|  
|-40442.1939107|52679.5599876|L00S\_1250W\_GT|  
|-40609.3157689|52495.8103203|L00S\_1500W\_GT|

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?

1. 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

gmt grdproject -Jm-81/46/1:1 -n+c -D100+ -R-81.6/-81.4/46.4/46.6 -F -C -I “tmp.nc?Conductivity[$index]” -V -Gmerge.grd

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.

1. 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.

Cheers

Max

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 28, 2020, 8:36pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/2 "2020-08-28T20:36:08Z")

</div>

Hi Max-

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.

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 28, 2020, 8:49pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/3 "2020-08-28T20:49:41Z")

</div>

> [@pwessel](#):
>
> -Jm-18/0/1:1

Should be -Jm-81/0/1:1 of course.

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [August 29, 2020, 11:18am UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/4 "2020-08-29T11:18:19Z")

</div>

Note that you can recover the 6.0 results by using proj4 syntax, _e.g._

```auto
echo -81.50638333 46.47973333 | gmt mapproject -J+proj=merc+lat_ts=46+lon_0=-81
-39226.1233521 4054500.74761

```

Maybe the `grdproject` works with it too. What does this do?

```auto
gmt grdproject -J+proj=merc+lat_ts=46+lon_0=-81 -n+c -D100+ -I “tmp.nc?Conductivity[$index]” -V -Gmerge.grd

```

---

<div class="post-metadata">

**Author:** ![pwessel](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/pwessel/32/30_2.png) [@pwessel](https://forum.generic-mapping-tools.org/u/pwessel)\
**Post date:** [August 30, 2020, 1:07am UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/5 "2020-08-30T01:07:06Z")

</div>

We fixed two issues in 6.1.0 that lead to this behavior. Master is now back to giving you your expected result; see [https://github.com/GenericMappingTools/gmt/pull/4120](https://github.com/GenericMappingTools/gmt/pull/4120) for details.

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.

---

<div class="post-metadata">

**Author:** ![mmoorkamp](https://avatars.discourse-cdn.com/v4/letter/m/f4b2a3/32.png) [@mmoorkamp](https://forum.generic-mapping-tools.org/u/mmoorkamp)\
**Post date:** [August 31, 2020, 12:36pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/6 "2020-08-31T12:36:48Z")

</div>

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]]

---

<div class="post-metadata">

**Author:** ![Joaquim](https://yyz1.discourse-cdn.com/flex047/user_avatar/forum.generic-mapping-tools.org/joaquim/32/16_2.png) [@Joaquim](https://forum.generic-mapping-tools.org/u/Joaquim)\
**Post date:** [August 31, 2020, 1:51pm UTC](https://forum.generic-mapping-tools.org/t/mapproject-grdproject-issues-and-changes-between-versions/823/7 "2020-08-31T13:51:29Z")

</div>

> [@mmoorkamp](#):
>
> grdproject [INFORMATION]: Converted to -J syntax = -Jm-81/0/1:1

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).