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

code123
  • 345
  • 4
  • 12

2 Answers2

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.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • 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 05 '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 05 '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 05 '16 at 18:12
  • target flips my map upside down. It is a map of Canada – code123 Feb 05 '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 05 '16 at 18:38
  • Here is my code: dat<- raster('2013_000_CaPA_continental_v2.3-analyse_final_10km_6h_APCP.nc')# read one layer only for illustration;

    print(dat) # print netcdf details; proj4string(dat) <- CRS"unknown";

    projj <- CRS("+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"); rx <- projectRaster(from=dat, crs=projj)

    cr <- crop(rx, extent(Prairie.Boundaries)) ; fr <- rasterize(shapfile, cr); Can <- mask(x=cr, mask=fr);;

     – code123 Feb 05 '16 at 18:40
  • I just replaced the larger file with a much smaller subset (just one timestep) from the rasterBrick. Hopefully you can use this to help me out – code123 Feb 05 '16 at 18:54
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.

Dmitry Baryshnikov
  • 3,221
  • 13
  • 17