Grdproject: convert X/Y grid to a lat/long grid

Q&A
1

Hello, I hope you can help. I am using gmt5.4.3:

I have a netcdf file which is in X/Y coordinate system: This is the information from grdinfo

input.nc: Command: xyz2grd shear_stress.xyz -I5000 -R-2000000/2500000/-1250000/3100000 -Ggrids/shear_stress_poly.nc
input.nc: Remark:
input.nc: Gridline node registration used
Grid file format: nf (# 18) GMT netCDF format (float) (COARDS-compliant) [DEFAULT]

  • x_min: -2000000 x_max: 2500000 x_inc: 5000 name: x nx: 901*
  • y_min: -1250000 y_max: 3100000 y_inc: 5000 name: y ny: 871*
  • z_min: 5000 z_max: 130000 name: z*
    scale_factor: 1 add_offset: 0

I want to convert to a latitude and longitude format: This is the information that was provided about the output original projection
project: albers
center_longitude=13.5
center_latitude=60.8
southern_standard_parallel=43
northern_standard_parallel=62
it has a scale of 5000m

I am using the command grdproject:

grdproject -Jl13.5/60.8/43/62/1:1 input.nc -I -C -F out.nc

However, the resultant latitude and longitude are not what I know there should be: i.e the input grid extends into the barents sea, ~ 75N but the final latitude and longitude limits are
Longitude: -12W and 100E
latitude: 46N to 68N

I am using the grdproject command correct given the information I was provided.

Thank!

Sarah

2

you are using -Jl, that is Lambert projection, but you says that it is a Albers projection, then you should use -Jb

3

ops good spot. Thanks. However, nearly the same issues:
the revised latitude and longitude in the out.nc are:

lat: 46.4N to 69N
lon: -13W to 108.8E

So still does not seem to relate to where I know the data is?

Sarah

4

using info from grdinfo:

cat << EOF > range.txt
-2000000 -1250000
2500000 -1250000
2500000 3100000
-2000000 3100000
EOF
gmt mapproject range.txt -Jb13.5/60.8/43/62/1:1 -Rd -I -C
0.712780291462 55.8100384291
29.3940832463 55.470337274
38.6651742095 70.4229781997
-6.91795688286 71.0017797023

Are these values close to what you expected?

Not sure how deal with the ‘it has a scale of 5000m’

using gmt 6.1.1

5

Hi Marcelo

I had not thought to try mapproject reading in the limits of the X/Y - I will give that a go.!

thank!

Sarah

6

Hi Marcelo

It is still not quite in the right place - it should include UK and Barents sea; I will double check the projection information I was provided.

Sarah

7

If you use a proj4 projection string in -J the -R option is not needed. In map|grdproject the -R extend may select either spherical or ellipsoidal approximations. With proj4 it’s always what those parameters set (default to ellipsoidal).

8

Hi Joaquim

I am a bit confused by what you are suggest: should I think remove the -R command from the grdproject or mapproject?

I have tried to find a way of defined proj4 in the man pages but nothing is listed.

I removed the -R command from mapproject commands: but I get an error to define -R:
mapproject range.txt -Jb13.5/60.8/43/62/1:1 -I -C

Thanks

9

The GMT docs explain on how to use proj4 strings but does not list them (they are thousands). You would need to consult the PROJ manuals to know how to write one from your projection parameters. Then with it you don’t need -R but with the -J projection syntax you need, which is a pain.

10

Hi Joaquim

I was just going by the grdproject information. I will go looking further into the GMT docs page.

Now that I know I am on the right track.

thanks for answering all my questions.

Sarah

11

$ gmt mapproject range.txt -Jaea/13.5/60.8/43/62/1:1 -I -C
mapproject [ERROR]: Must specify -R option

Appears that mapproject with proj need -R
using gmt 6.2.0_50eafd2_2021.02.10

12

You are not using -J+proj=… That’s using mappoject with proj

13

Hi

Thanks for the suggestions. But as I only have version 5.4.3 not version6. I am working on HPC cluster and I do not have the access to GMT6.

From GMT6 docs; I can see all the many proj codes: -Jproj|EPSG.
I have tried to use the listed codes for GMT5:

grdproject -Jaea/13.5/60.8/43/62/1:1 input.nc -I -C -F -Gout_latlon.nc

however, I get an error message
grdproject: Internal Error = GMT_MAP_BAD_LAT_MAXgrdproject (GMT_grdproject): North is outside -90 to +90 degree range

I will see if I can upgrade to gmt6; but I might have to look into a different software.

Sarah

14

You could consider doing the pre and post processing in more simple computer than the HPC.

The following commands produce the same results:

gmt mapproject range.txt -Jb13.5/60.8/43/62/1:1 -Rd -I -C -Fe

proj=‘-J+proj=aea+lon_0=13.5+lat_0=60.8+lat_1=43+lat_2=62+ellps=WGS84+units=m’
gmt mapproject range.txt $proj -I

Just to you have some more things to try.

15

I would love to be install things on a local machine - but I do not have that option with this newer version.
Thanks for so much for trying the non-GMT6 version I will keep trying with these extra commands this is the output I get; - which i am going for is the four corners of a lat/lon square.

-13.2872884533 46.4046346244
46.1998312515 44.6351070531
103.977549608 73.0906845128
-71.4534740956 78.0367484823

16

If you have Anaconda or miniconda installed on your HPC, you can install GMT6 using the following command:

conda install -c conda-forge gmt

Then GMT6 should be installed in your own home directory.

17

Hi, never knew could try that. great tip!
giving it a go now and so far - not crashed!

Sarah