How to project x,y coordinates to lat/lon in netcdf file

Question

I have downloaded the velocity field of the Greenland ice sheet from the CCI website as a NetCDF file. However, the projection is given as (see below, where x ranges between [-639750,855750] and y [-655750,-3355750])

How can I project these data to actual lat/lon coordinates in the NetCDF file? Thanks already! For the ones interested: the file can be downloaded here: http://products.esa-icesheets-cci.org/products/downloadlist/IV/

Variables:
    crs                                
           Size:       1x1
           Dimensions: 
           Datatype:   int32
           Attributes:
                       grid_mapping_name                     = 'polar_stereographic'
                       standard_parallel                     = 70
                       straight_vertical_longitude_from_pole = -45
                       false_easting                         = 0
                       false_northing                        = 0
                       unit                                  = 'meter'
                       latitude_of_projection_origin         = 90
                       spatial_ref                           = 'PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",GEOGCS["WGS 84",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","3413"]]'
    y                                  
           Size:       5401x1
           Dimensions: y
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'Y'
                       long_name     = 'y coordinate of projection'
                       standard_name = 'projection_y_coordinate'
    x                                  
           Size:       2992x1
           Dimensions: x
           Datatype:   double
           Attributes:
                       units         = 'm'
                       axis          = 'X'
                       long_name     = 'x coordinate of projection'
                       standard_name = 'projection_x_coordinate'

what has Python got to do with this? – WiseDev Jul 02 '19 at 12:03 — WiseDev, Jul 02 '19 at 12:03

score 6 · Accepted Answer · answered Jul 05 '19 at 05:29

If you want to transform the whole grid from its native Polar Stereographic coordinates to a geographic (longitude by latitude) grid, you'll probably want to use a tool like gdalwarp. I don't think that's the question you're asking, though.

If I'm reading your question correctly, you want to pick points out of the file and locate them as lon/lat coordinate pairs. I'm assuming that you know how to get a location as an XY pair out of your netCDF file, along with the velocity values at that location. I'm also assuming that you're doing this in Python, since you put that tag on this question.

Once you've got an XY pair, you just need a function (with a bunch of parameters) to transform it to lon/lat. You can find that function in the pyproj module.

Pyproj wraps the proj4 C library, which is very widely used for coordinate system transformations. If you have an XY pair in projected coordinates and you know the definition of the projected coordinate system, you can use pyproj's transform function like this:

import pyproj

# Output coordinates are in WGS 84 longitude and latitude
projOut = pyproj.Proj(init='epsg:4326')

# Input coordinates are in meters on the Polar Stereographic 
# projection given in the netCDF file
projIn = pyproj.Proj('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 
    +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs ',
    preserve_units=True)

# here is a coordinate pair near the middle of your data set
x, y = 0.0, -2000000

# transform x,y to lon/lat
lon, lat = pyproj.transform(projIn, projOut, x, y)

# answer: lon = -45.0; lat = 71.6886

... and there you go. Note that the output longitude is -45.0, which should give you a nice warm feeling, since the input X coordinate was 0, and -45.0 is the central meridian of the data set's projection. If you want your answer in radians instead of degrees, set the radians kwarg in the transform function to True.

Now for the hard part, which is actually the thing you do first -- defining the projIn and projOut that are used as arguments for the transform function. These are in the input and output coordinate systems for the transformation. These are Proj objects, and they hold a mess of parameters for the coordinate system transformation equations. The proj4 developers have encapsulated them all in a tidy set of functions and the pyproj developers have put a nice python wrapper around them, so you and I don't have to keep track of all the details. I will be grateful to them for all the days that remain to me.

The output coordinate system is trivial

projOut = pyproj.Proj(init='epsg:4326')

The pyproj library can build a Proj object from an EPSG code. 4326 is the EPSG code for WGS 84 lon/lat. Done.

Setting projIn is harder, because your netCDF file defines its coordinate system with a WKT string, which (I'm pretty sure) can't be read directly by proj4 or pyproj. However, pyproj.Proj() will take a proj4 parameter string as an argument. I've already given you the one you need for this operation, so you can just take my for for it that this

+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

is the equivalent of this (which is copied directly from your netCDF file):

PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
  GEOGCS["WGS 84",
    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"]],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Polar_Stereographic"],
  PARAMETER["latitude_of_origin",70],
  PARAMETER["central_meridian",-45],
  PARAMETER["scale_factor",1],
  PARAMETER["false_easting",0],
  PARAMETER["false_northing",0],
  UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
AUTHORITY["EPSG","3413"]]'

If you want to be able to do this more generally, you'll need another module to convert WKT coordinate system definitions to proj4 parameter strings. One such module is osgeo.osr and there's an example program at this blog post that shows you how to do that conversion.

Hi Tom! Thanks, this really helped a lot! So for future cases, if I need to to do whole projection, would this command do the job with gdalwarp? gdalwarp -overwrite -s_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs EPSG:4326 greenland_polarstereo.nc greenland_latlon.nc — Yoni Verhaegen, Jul 05 '19 at 08:46
You're on the right track. It's been a while since I used gdalwarp enough to remember its many options, though. I suggest playing with a smaller file to figure out what the options do. You might not need the -s_srs option, for example, because gdalwarp may read the input CRS from the netCDF file. IMPORTANT!=>Specify an output resolution with the -tr option, and resampling method with the -r option. You can't change the CRS of a grid without resampling, and that means you're interpolating between the original data values. Choose the interpolation resolution and method purposefully. — TomEvans, Jul 06 '19 at 03:45

How to project x,y coordinates to lat/lon in netcdf file

1 Answers1