Plot Polar Stereographic Grid

Hi, all

I have a problem to plot the Netcdf data in Sea Ice Polar Stereographic North projection which is used in present NSIDC dataset.

The Sea Ice Polar Stereographic North information is here: WGS 84 / NSIDC Sea Ice Polar Stereographic North - EPSG:3413

I tried to convert the XY coordinate to geo-coordinate, but seems not correct. And the data values are all NAN. I am not sure why?

Here is the pyGMT code:

!gmt grdproject -JEPSG:3413 -R0/360/50/90 -C -I

fig = pygmt.Figure()

fig.basemap(region="0/360/40/90", projection="s0/90/71/1:60000000", frame="a30f")
pygmt.makecpt(cmap="haxby", series=[-1, 1])
# grid=''
fig.colorbar(frame=["af", "y+lm"])

And the original NC file: (918.3 KB)

The data could be plotted by Panoplay, which is correct:



Have not checked or tested, but is the dataset in geodetic or projected coordinates? If latter, then you just plot it as cartesian data;




This is because they have already been projected (ref. Epsg:…).

The NSIDC files are already projected to the Sea Ice Polar Stereographic X-Y coordinate.
Here is the info of file

!gmt grd2xyz | head

grd2xyz [WARNING]: "y", NetCDF: Variable not found
	If something bad happens later, try importing via GDAL.
grd2xyz [WARNING]: "x", NetCDF: Variable not found
	If something bad happens later, try importing via GDAL.
0.498355263158	447	-9999
1.49506578947	447	-9999
2.49177631579	447	-9999
3.48848684211	447	-9999
4.48519736842	447	-9999
5.48190789474	447	-9999
6.47861842105	447	-9999
7.47532894737	447	-9999
8.47203947368	447	-9999
9.46875	447	-9999

This file has some -9999 values, but in the center the values are valid.

The was in geo coordinate using the grdproject. So I use the projection="s0/90/71/1:60000000" to plot the file.

I find this in the grdproject document:

Nodes not constrained by input data are set to NaN.

Is this the reason for all NaN? What should I do to avoid the valid data to set to NaN?

Not that I understand a lot, but your grid X and Y seems scaled. See more in the code below.
I managed to plot it with the following shell script:

    # flip the grid vertically
    gmt grdedit "" -Ev
    # set ice thickness below 0 to NaN
    gmt grdclip -Sb0/NaN 
    # plot
    gmt begin XXX png 
      gmt grdimage -Jx1:6000 -Xc -Yc
      # coastline, stereographic projection starting at -45 degrees
      gmt coast -R-180/180/40/90 -Js-45/90/40/1:100000000 -Bafg15 -Dc -W0.25p -X -Y
      gmt colorbar -DJRM+w6.5c/0.5c+o1c/0+mc -F+p+i -Bxa -By+lm
    gmt end

the result:

Note that your netCDF file is not COARDS compliant so you probably need to read it via GDAL. GMT expects all grids and cubes to be following the COARDS convention.

Thanks very much. You gave a good solution. But at present, I can not produce the same figure as yours.
If I keep the -Xc -Yc in grdimage sentence, the figure is very narrow and long and no data were plotted. If I remove the -Xc -Yc, I can see the data with wrong location.

Maybe its related to GMT version (GMT 6.1.0)?

I will give more try.
Thanks again.

My first suggestion was wrong.
I assumed the grid center is at the pole. This was wrong, it is not. This is why northern coast of Greenland and Canadian arctic islands did not match the southern edge of ice on the grid.

The RDEFT4_2010102 grid contains datasets called “lon” and “lat”, as well as the ice variables. the X and Y are scaled, but it is not clear to me how to place that grid on a bigger map. The metadata (text abstract) says it is a 25-km grid, but the Cartesian coordinates are relative to one of the corners of the grid which is set to zero.

Extracting the sea_ice_thickness “subdataset” with GDAL made no difference. The result was the same Cartesian grid as produced just with grdedit.

another approach gave a better result. I extracted (manually) low left and top right corner coordinates of the grid using GDAL and converted it to text with xyz2grd fom subdatasets “lon” and “lat” and used these to produce a stereographic projection map with a rectangular ROI with the map size matching the grid size.


rm gmt*
# flip the grid vertically
gmt grdedit "" -Ev
# set ice thickness below 0 to NaN
gmt grdclip -Sb0/NaN 
# plot
WIDTH=10.0 # 10 cm
# 447 is grid height, 303 is grid width; -JX needs height indicated separately if height does not match width
HEIGHT=$(awk -v WIDTH="$WIDTH" 'BEGIN {print 447/303*WIDTH}')
gmt begin XXX png 
  gmt grdimage -JX${WIDTH}c/${HEIGHT}c
  # coastline, stereographic projection starting at -45 degrees
  gmt coast -R279.28500/34.05146/102.37037/31.487499r -JS-45/90/${WIDTH}c -Bxa90-45g15 -Byg15 -Di -W0.25p #-Gwhite
  gmt colorbar -DJRM+w5c/0.5c+o1.5c/0+mc -F+p+i -Bxa -By+lm
gmt end

the result:

Great! I can re-produce the figure. Thank you.

Let’s check the data structure.

ncdump | head -200
netcdf RDEFT4_20101021 {
	x = 304 ;
	y = 448 ;
	float sea_ice_thickness(y, x) ;
		sea_ice_thickness:units = "Meters" ;
		sea_ice_thickness:long_name = "Sea ice thickness" ;
	float snow_depth(y, x) ;
		snow_depth:units = "Meters" ;
		snow_depth:long_name = "Snow depth" ;
	float snow_density(y, x) ;
		snow_density:units = "Kg / Meters^3" ;
		snow_density:long_name = "Snow density" ;
	float lat(y, x) ;
		lat:units = "Degrees" ;
		lat:long_name = "Latitude" ;
	float lon(y, x) ;
		lon:units = "Degrees" ;
		lon:long_name = "Longitude" ;
	float freeboard(y, x) ;
		freeboard:units = "Meters" ;
		freeboard:long_name = "Ice freeboard" ;
	float roughness(y, x) ;
		roughness:units = "Meters" ;
		roughness:long_name = "Ice surface roughness" ;
	float ice_con(y, x) ;
		ice_con:units = "Percent x 100" ;
		ice_con:long_name = "Sea ice concentration" ;

// global attributes:
		:Title = "CryoSat-2 sea ice thickness and ancillary data" ;
		:Abstract = "This data set contains monthly averaged Arctic sea ice thickness estimates and ancillary data. The primary data set used in the production of these data come from the ESA CryoSat-2 satellite. Sea ice freeboard is determined from CryoSat-2 using the method described in the Reference section below. In brief, this method uses a physical model to determine the best fit to each CryoSat-2 waveform.The fitted waveform is used to determine the retracking correction and also allows determination of the surface roughness within the footprint. For sea ice floes, the dominant backscattering layer is taken to be from the sea ice surface and thus sea ice freeboard is here defined as the height of the ice layer above the local sea surface. The DTU10 MSS is subtracted from each elevation measurement and the elevations from leads and sea ice floes are placed onto a 25 km polar stereographic grid. Sea ice freeboard is then determined by subtracting the gridded sea surface elevation from the gridded sea ice floe elevation. Snow depth is constructed from a modified Warren climatology of snow depth on sea ice. Sea ice thickness is retrieved assuming hydrostatic balance and nominal densities of snow, ice, and water. Retrievals are only done when the sea ice concentration is at least 70%. Sea ice concentration is from the near real time DMSP SSMI_S daily polar gridded data set with the pole hole set to a constant value of 100%." ;
		:Projection = "CryoSat-2 elevation data have a nominal footprint size of 380 m x 1,650 m and are referenced to the WGS-84 ellipsoid. The derived CryoSat-2 freeboard data, and all ancillary data, have been gridded to the 25 km polar stereographic SSM/I grid. The center latitude and longitude for each grid point are provided with the data." ;

We find that the x and y values are just the index of the two dimensions. They are not the projected coordinates in the Polar Stereographic projection and their unit is not meter (just index).

The lat and lon are stored in the variables as with index of (y,x). As we dump the variables using grd2xyz, the true value of the latitude and longitude could be seen. However, the geo-coordinates not in regular grid nodes.

The GMT can not plot this kind of data directly because of the special data structure. GMT can not read the true coordinate no matter the geo or xy forms.

After a coordinates transform as adopted by mkononets, the data can be plotted.

So, maybe the GMT developers could consider adding a new feature about this. It will be very useful to plot this projection which is widely used in polar area.

gdal_transform can be used to transform the grid X,Y indices into meters using on the projection info from A Guide to NSIDC's Polar Stereographic Projection | National Snow and Ice Data Center

the end result is similar to the example 28, except the linear scale annotation on the map axes make no sense when using stereographic projection.


rm -f gmt.*

gdal_translate -of NETCDF -a_ullr -3850000 -5350000 3750000 5850000 -a_srs EPSG:3413 NETCDF:"":sea_ice_thickness
# set ice thickness below 0 to NaN
gmt grdclip -Sb0/NaN 
# define map scale
gmt begin RDEFT4_20101021_sea_ice_thickness1 png 
  # plot Cartesian grid data, in meters, scale as defined above
  gmt grdimage -Jx$SCALE
  # coastline, stereographic projection, origin at -45 degrees lon/ 90 degrees lat. 
  # Scale (defined above) valid at a standard latitude 70 degrees (no distortion at 70 degreed north)
  # -Js-45/90/70/1:75000000
  gmt coast -Js-45/90/70/$SCALE -Bxa90-45g15 -Byg15 -Di -W0.25p #-Gwhite
  gmt colorbar -DJRM+w5c/0.5c+o1.5c/0+mc -F+p+i -Bxa -By+lm
gmt end

@mkononets Thanks a lot for providing another excellent method! It is very helpful to me.

I read that gdal has Python bindings too (import gdal; you might need to install it first)
If you manage to put everything together in a single Python script, calling both gdal and GMT from Python, publish it here.

@mkononets, I just put your shell code to the jupyter notebook. See:

The gdal and GMT executable programs not wrapped by pygmt could be called in Jupyter using ! before programs. Or, you can use the python function os.system in pure python script like:

os.system(" gdal_translate -of NETCDF -a_ullr -3850000 -5350000 3750000 5850000 -a_srs EPSG:3413 NETCDF:'':sea_ice_thickness ")

This is useful for data plotting and processing, but not so pure python because we use some programs by system calling.

Sure, we can translate them to pure python script by importing gdal lib. But, I think the jupyter is better for pygmt and the system calling is acceptable.

Yes running gdal and gmt executables from a Jupiter notebook is acceptable for just solving this task. But there is nothing interesting to me about this approach. Jupiter is not even necessary for performing such task. It is much more interesting to get the data into Python taking into account the vast possibilities of transforming and processing data by means of Python (gdal, pyGMT, and what not).

Just thinking loudly.

for unclear (to me) reasons the previous version of grid processing commands (gdal_translate + gmt grdclip) resulted in wrong Y scaling of the resulting grid, y_inc: 24888.39 after running gdal_translate and y_inc: 24777.28 after gmt grdclip

an updated variant of gdal_translate command (below) replaces those two, preserves vertical grid y_inc: 25000.00 and assigns the correct NoData value.

gdal_translate --config GDAL_NETCDF_BOTTOMUP NO \
  -of netCDF -a_ullr -3850000 5850000 3750000 -5350000 -a_srs EPSG:3413 -a_nodata -9999 \
  NETCDF:"":sea_ice_thickness \

it changes grid registration from gridline to pixel, but seems to preserve both the correct scale and the boundaries.