Struggling with converting Cartesian X,Y,Z coordinates to Lat,Lon,Elevation with map project

Q&A
1

Hi there,
I’ve got a netCDF file including bed elevation for the Antarctic region, with the data in cartesian coordinates. I can extract the x, y and bed topography data from the netCDF, for instance, so as to build a 3-column text file to give to mapproject.
I just cannot figure out how to use this properly; have tried many options but do not get any sensible output.
Typically like: gmt mapproject XYZ.txt -Js0/-90/2.8i/-30 -R-180/180/-90/-30 -I -E > LATLONZ.txt
Example of the first 20 lines of the XYZ file below, as well as the header information in the netCDF file.
Anyone with some gook geodesy skills to help?
Thanks,
David

#######
-3333000.00 3333000.00 -5915.61719
-3333000.00 3328000.00 -5962.10742
-3333000.00 3323000.00 -5959.86377
-3333000.00 3318000.00 -5866.03711
-3333000.00 3313000.00 -5661.64453
-3333000.00 3308000.00 -5622.04150
-3333000.00 3303000.00 -5500.58740
-3333000.00 3298000.00 -5393.52588
-3333000.00 3293000.00 -5370.03662
-3333000.00 3288000.00 -5484.46045
-3333000.00 3283000.00 -5668.95605
-3333000.00 3278000.00 -5865.19629
-3333000.00 3273000.00 -5878.93457
-3333000.00 3268000.00 -5886.30225
-3333000.00 3263000.00 -5934.98633
-3333000.00 3258000.00 -5999.31201
-3333000.00 3253000.00 -6071.00391
-3333000.00 3248000.00 -6139.77783
-3333000.00 3243000.00 -6057.91699
-3333000.00 3238000.00 -6030.22168
etc…

####### header of the netCDF file.
netcdf BedMachineAntarctica_2020-07-15_v02 {
dimensions:
x = 13333 ;
y = 13333 ;
variables:
char mapping ;
mapping:grid_mapping_name = “polar_stereographic” ;
mapping:latitude_of_projection_origin = -90. ;
mapping:standard_parallel = -71. ;
mapping:straight_vertical_longitude_from_pole = 0. ;
mapping:semi_major_axis = 6378273. ;
mapping:inverse_flattening = 298.27940504282 ;
mapping:false_easting = 0. ;
mapping:false_northing = 0. ;
int x(x) ;
x:long_name = “Cartesian x-coordinate” ;
x:standard_name = “projection_x_coordinate” ;
x:units = “meter” ;
int y(y) ;
y:long_name = “Cartesian y-coordinate” ;
y:standard_name = “projection_y_coordinate” ;
y:units = “meter” ;
float bed(y, x) ;
bed:long_name = “bed topography” ;
bed:standard_name = “bedrock_altitude” ;
bed:units = “meters” ;
bed:grid_mapping = “mapping” ;
bed:source = “IBCSO and Mathieu Morlighem” ;
bed:_FillValue = -9999.f ;
short geoid(y, x) ;
geoid:long_name = “EIGEN-EC4 Geoid - WGS84 Ellipsoid difference” ;
geoid:standard_name = “geoid_height_above_reference_ellipsoid” ;
geoid:units = “meters” ;
geoid:grid_mapping = “mapping” ;
geoid:geoid = “eigen-6c4 (Forste et al 2014)” ;

// global attributes:
:Conventions = “CF-1.7” ;
:Title = “BedMachine Antarctica” ;
:Author = “Mathieu Morlighem” ;
:version = “15-Jul-2020 (v2.0)” ;
:nx = 13333. ;
:ny = 13333. ;
:Projection = “Polar Stereographic South (71S,0E)” ;
:proj4 = “+init=epsg:3031” ;
:sea_water_density\ (kg\ m-3) = 1027. ;
:ice_density\ (kg\ m-3) = 917. ;
:xmin = -3333000 ;
:ymax = 3333000 ;
:spacing = 500 ;
:no_data = -9999. ;
:license = “No restrictions on access or use” ;
:Data_citation = “Morlighem M. et al., (2019), Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nature Geoscience (accepted)” ;
:Notes = “Data processed at the Department of Earth System Science, University of California, Irvine” ;
}
##############

2

Try using the EPSG code (3031 for Antarctic Polar Stereographic). The below seems to give sensible-ish results, I’m not sure why the mapproject -E option doesn’t work as expected with the inverse -I option.

gmt mapproject XYZ.txt -JEPSG:3031 -R-180/180/-90/-60 -I > LONLATZ.txt

produces the following (note that it’s Longitude, Latitude, Elevation). Not sure if the Elevation is correct, but at least it looks better than before…

-45	-48.4643831064	-5915.61719
-45.0430083837	-48.4929562939	-5962.10742
-45.0860813347	-48.5215132292	-5959.86377
-45.1292189499	-48.5500538498	-5866.03711
-45.1724213258	-48.5785780935	-5661.64453
-45.2156885592	-48.6070858975	-5622.0415
-45.2590207468	-48.6355771991	-5500.5874
-45.3024179853	-48.6640519356	-5393.52588
-45.3458803713	-48.6925100438	-5370.03662
-45.3894080017	-48.7209514607	-5484.46045
-45.4330009729	-48.7493761231	-5668.95605
-45.4766593818	-48.7777839675	-5865.19629
-45.5203833248	-48.8061749305	-5878.93457
-45.5641728987	-48.8345489484	-5886.30225
-45.6080282001	-48.8629059574	-5934.98633
-45.6519493255	-48.8912458938	-5999.31201
-45.6959363715	-48.9195686934	-6071.00391
-45.7399894347	-48.9478742921	-6139.77783
-45.7841086116	-48.9761626256	-6057.91699
-45.8282939987	-49.0044336294	-6030.22168

And just to busybody a bit, is there a reason to convert the BedMachine Antarctica NetCDF grid to an XYZ file before reprojecting? Just wondering if using grdproject might be better, though I’m not sure what you’re trying to do exactly.

3

Oh, well, that works actually!
I had first to install GMT6 (was still on 5) so that the -JEPSG:3031 was recognised.
I do not work directly from the netCDF file because it contains multiple variables and I’m not sure how to tell to grdproject which one he has to read.
Thanks a lot for your help.
Have a nice WE

4

Yay, glad that it worked :smile_cat:

FYI, you can use file.nc?varname to access a single NetCDF variable as per 3. General Features — GMT 6.2.0 documentation. E.g. gmt grdinfo BedMachineAntarctica_2020-07-15_v02.nc?bed will give:

BedMachineAntarctica_2020-07-15_v02.nc: Title: bed topography
BedMachineAntarctica_2020-07-15_v02.nc: Command: 
BedMachineAntarctica_2020-07-15_v02.nc: Remark: 
BedMachineAntarctica_2020-07-15_v02.nc: Gridline node registration used [Cartesian grid]
BedMachineAntarctica_2020-07-15_v02.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
BedMachineAntarctica_2020-07-15_v02.nc: x_min: -3333000 x_max: 3333000 x_inc: 500 name: Cartesian x-coordinate [meter] n_columns: 13333
BedMachineAntarctica_2020-07-15_v02.nc: y_min: -3333000 y_max: 3333000 y_inc: 500 name: Cartesian y-coordinate [meter] n_rows: 13333
BedMachineAntarctica_2020-07-15_v02.nc: v_min: 0 v_max: 0 name: bed topography [meters]
BedMachineAntarctica_2020-07-15_v02.nc: scale_factor: 1 add_offset: 0
BedMachineAntarctica_2020-07-15_v02.nc: format: netCDF-4 chunk_size: 953,953 shuffle: on deflation_level: 1

I tried running gmt grdproject on BedMachine Antarctica but it took a long time and didn’t finish. If you have a smaller study region (e.g. a single glacier/ice stream), it might work though. Or maybe just stick with mapproject :upside_down_face:

5

Thanks for the tip on how to select one variable from a netCDF file.
Looks like I’m not very good at digging into the documentation (!)
Regards,
David