2

I have a netCDF file and would like to convert it to LAT-LON coordinates. 1) I do not know the projection in the file but from the website of the data provider, the data is on a polar-stereographic (PS) grid covering North America and adjacent waters with a 10 km resolution at 60 degrees north. 2) How can I project to latlon using a suitable coordinate system for the north pole?

Here is my sample file:

https://www.dropbox.com/s/lgni1s5tgu1sbwh/data.nc?dl=0

  • Any chance of something smaller than 750Mb to download? Maybe just a few coordinate points in the given coordinate system? – Spacedman Feb 5 '16 at 16:35
  • @Spacedman I just added a smaller file size `3MB – code123 Feb 5 '16 at 18:54
  • Unfortunately, GDAL 2.1.0 can not read the reduced file. – AndreJ Aug 16 '16 at 9:07
1

The projection looks like it matches this custom entry in http://spatialreference.org:

PROJCS["Stereographic_North_Pole",
    GEOGCS["GCS_Coordinate System imported from GRIB file",
        DATUM["D_unknown",
            SPHEROID["Sphere",6371229,0]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Stereographic_North_Pole"],
    PARAMETER["standard_parallel_1",60],
    PARAMETER["central_meridian",249],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

http://spatialreference.org/ref/sr-org/8237/

Depending on the software you have at hand, that should be usable in one format or another.

If you want to transform a shapefile in lat-long to this grid system, then try this:

require(rgdal)
target = "+proj=stere +lat_0=90 +lat_ts=60 +lon_0=249 +k=1 +x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=m +no_defs"
shgrid = spTransform(shlatlong, CRS(target))

plot the shgrid object to make sure. You might get problems if the shapefile crosses the dateline or rounds the pole...

Not knowing how you've loaded your netCDF into R, I can't be sure how to crop it, but if its a raster then mask from the raster package probably does the trick.

  • it is true that It might not be a good idea to transform a grid over the pole to lat-long but I have a shapefile in latlong which I would like to use in cropping the netcdf file. Here are the bbox coordinates of the shapefile: xmin -141.01807 xmax -52.61941 ymin 41.68144 ymax 83.13550 Can you show me how to use these coordinates to crop the file attached here? I can also include the shapefile in the attachment if you need it. – code123 Feb 5 '16 at 17:41
  • Do you have QGIS, ArcGIS, R, something else, or do you want to do all this with GDAL command-line tools? – Spacedman Feb 5 '16 at 18:03
  • I have R. I can use raster and rasterVis packages for doing it if you give me directions. gdaltransform -s_srs SR-ORG:8237 -t_srs epsg:5937 2013_000_CaPA_continental_v2.3-analyse_final_10km_6h_APCP.nc 2013.tif gdal does not recognize SR-ORG:8237. It issues an error. – code123 Feb 5 '16 at 18:12
  • target flips my map upside down. It is a map of Canada – code123 Feb 5 '16 at 18:33
  • Ooh. Possibly because the grid system reference has Y increasing southwards (on the Canadian side of the pole anyway). Sorry, but I can't download a 750Mb file to play with this. – Spacedman Feb 5 '16 at 18:38
0

To reproject raster into another spatial reference the gdalwarp utility can be used. The command can be looks like:

gdalwarp -t_srs 'EPSG:4326' -overwrite \
NETCDF:"your file.nc":your_subdataset subdataset_4326.tif

Choose another projection if you wish.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.