Using Arctic DEM

andreiCCI · January 14, 2020, 1:55am

Hello,
I am trying to plot Arc DEM using GMT 6.0 installed using conda on Mac OS.
I downloaded 1 kilometer ArcticDEM Mosaic arcticdem_mosaic_1km_v3.0.tif file DEM. According to the web site Mosaicked DEMs are provided in 32-bit GeoTIFF format.

1. This is file info:

gmt grdinfo arcticdem_mosaic_1km_v3.0.tif
arcticdem_mosaic_1km_v3.0.tif: Title: Grid imported via GDAL
arcticdem_mosaic_1km_v3.0.tif: Command:
arcticdem_mosaic_1km_v3.0.tif: Remark:
arcticdem_mosaic_1km_v3.0.tif: Pixel node registration used [Cartesian grid]
arcticdem_mosaic_1km_v3.0.tif: Grid file format: gd = Import/export through GDAL
arcticdem_mosaic_1km_v3.0.tif: x_min: -4000000 x_max: 3400000 x_inc: 1000 name: n_columns: 7400
arcticdem_mosaic_1km_v3.0.tif: y_min: -3400000 y_max: 4100000 y_inc: 1000 name: n_rows: 7500
arcticdem_mosaic_1km_v3.0.tif: z_min: -718.904541016 z_max: 5912.12158203 name: z
arcticdem_mosaic_1km_v3.0.tif: scale_factor: 1 add_offset: 0
+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

I tried to convert it to GMT native netcdf format

gmt grdconvert arcticdem_mosaic_1km_v3.0.tif -Garcticdem_mosaic_1km_v3.0.nc -V

This is new Netcdf file info

gmt grdinfo arcticdem_mosaic_1km_v3.0.nc
arcticdem_mosaic_1km_v3.0.nc: Title: Produced by grdconvert
arcticdem_mosaic_1km_v3.0.nc: Command: grdconvert arcticdem_mosaic_1km_v3.0.tif -Garcticdem_mosaic_1km_v3.0.nc -V
arcticdem_mosaic_1km_v3.0.nc: Remark:
arcticdem_mosaic_1km_v3.0.nc: Pixel node registration used [Cartesian grid]
arcticdem_mosaic_1km_v3.0.nc: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7
arcticdem_mosaic_1km_v3.0.nc: x_min: -4000000 x_max: 3400000 x_inc: 1000 name: n_columns: 7400
arcticdem_mosaic_1km_v3.0.nc: y_min: -3400000 y_max: 4100000 y_inc: 1000 name: n_rows: 7500
arcticdem_mosaic_1km_v3.0.nc: z_min: -718.904541016 z_max: 5912.12158203 name: z
arcticdem_mosaic_1km_v3.0.nc: scale_factor: 1 add_offset: 0
arcticdem_mosaic_1km_v3.0.nc: format: netCDF-4 chunk_size: 130,130 shuffle: on deflation_level: 3
PROJCS[“unknown”,
GEOGCS[“unknown”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.257223563,
AUTHORITY[“EPSG”,“7030”]],
AUTHORITY[“EPSG”,“6326”]],
PRIMEM[“Greenwich”,0,
AUTHORITY[“EPSG”,“8901”]],
UNIT[“degree”,0.0174532925199433,
AUTHORITY[“EPSG”,“9122”]]],
PROJECTION[“Polar_Stereographic”],
PARAMETER[“latitude_of_origin”,70],
PARAMETER[“central_meridian”,-45],
PARAMETER[“false_easting”,0],
PARAMETER[“false_northing”,0],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,“9001”]],
AXIS[“Easting”,SOUTH],
AXIS[“Northing”,SOUTH]]

I tried to make some changes:

gmt grdedit arcticdem_mosaic_1km_v3.0.nc -R-180/180/60/90 -D+t"Arctic DEM Mosaic 500 m" -T -fg

and made a figure
gmt grdimage -Js0/90/12c/60 -R-180/180/60/90 arcticdem_mosaic_1km_v3.0.nc -pdf arcticdem1km
The grid is distorted.

Any help is greatly appreciated.
Thank you,
Andrei

pwessel · January 14, 2020, 2:09am

The grid is a cartesian rectangle in projected meters. grdedit does not undo any projection. It only replaced -R-4000000/3400000/-3400000/3400000 with -R-180/180/60/90. You need to run grdproject -I to undo the projection you really want to obtain a lon-lat grid. See the grdproject documentation. Alternatively, if you are OK with making the map in the stereographic projection you can plot it using a cartesian scale and overlay a geographic basemap using the same projection.

andreiCCI · January 14, 2020, 2:21pm

Dear Paul,
Thank you!

Would you be able to provide a full set of commands to convert ArcDEM GeoTiF file to the same format as @earth_relief_02m grid with the same names for the variables?
For me it would be easier to follow grdproject documentation with a working example.
I am trying to use Arctic DEM instead of @earth_relief_02m grd to generate elevation map for Greenland for the following figure below.
Very best,
Andrei

Joaquim · January 14, 2020, 5:17pm

Don’t know if it’s really a good idea to convert to geogs since it cover a pole but it should be something like

gmt grdproject arcticdem_mosaic_1km_v3.0.nc -I -Garctidem_geog.nc

andreiCCI · January 14, 2020, 6:08pm

Dear Joaquim,
Thank you.
I tried it and got the following message:
gmt grdproject arcticdem_mosaic_1km_v3.0.nc -I -Garctidem_geog.nc
grdproject [ERROR]: Syntax error: Must specify a map projection with the -J option

Joaquim · January 14, 2020, 6:29pm

What GMT version? With GMT6 it should not need a -J because your file is already georeferenced (I think).

andreiCCI · January 14, 2020, 6:42pm

Dear Joaquim,
It is gmt 6 GMT - The Generic Mapping Tools, Version 6.0.0 [64-bit] [4 cores]
installed using conda on Mac OS X

Joaquim · January 15, 2020, 3:43pm

I’m not able to reverse project that grid with nothing. This is the hard GMT way but …

grdproject arcticdem_mosaic_1km_v3.0.nc -I -Glixo.nc -Js-45/90/70/1:1 -V -C -R-180/180/0/90

even gdalwarp misses it

gdalwarp -t_srs "+proj=lonlat" arcticdem_mosaic_1km_v3.0.nc lixo.nc

andreiCCI · January 15, 2020, 4:19pm

Dear Joaquim,
Thank you again.
I clearly have no chance …
I found this data-description info online
https://www.pgc.umn.edu/guides/arcticdem/data-description/
Very best,
Andrei