Search code examples
rrasternetcdf

Convert UTM to latlon rasterBrick R

library(raster)
library(ncdf)
library(rgdal)

I am facing some issues converting UTM to latlon in a rasterBrick object. My code is below:

Example file can be found here: https://www.dropbox.com/s/yv4opmrx1v4whpt/2013_000_CaPA_continental_v2.3-analyse_final_10km_6h_APCP.nc.7z?dl=0

rasnc<- brick('file.nc', varname = "APCP_surface")

rasnc
#class       : RasterBrick 
#dimensions  : 824, 935, 770440, 1460  (nrow, ncol, ncell, nlayers)
#resolution  : 10000, 10000  (x, y)
#extent      : -5000, 9345000, -5000, 8235000  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 

wgs84.p4s <- "+proj=longlat +datum=NAD83 +ellps=GRS80 +no_defs"
# reprojecting
rx <- projectRaster(from=rasnc, crs=wgs84.p4s,method="ngb")

The above does not convert from UTM to latlon. I even tried to write the file to GeoTIFF and then re-project but it still did not work.

writeRaster(rasnc, file="myfil.tif", format="GTiff", overwrite=TRUE)

Basically I am trying to:

  • read a netCDF file.nc using brick

  • convert the UTM coordinates to latlon

  • crop the RasterBrick object to ext <- extent(-141.01807,-52.61941,41.68144.0,83.13550)

  • display any of the layers using levelplot function.

The time series is 6-hourly from Jan-01-2013 to dec-01-2013

It is true my use of the older raster package version was the source of the problem as suggested by @RobertH. rasnc has no coord.ref. From the webpage (http://weather.gc.ca/grib/grib2_RDPA_ps10km_e.html) of the data modellers, it is stated a polar-stereographic (PS) grid covering North America and adjacent waters with a 10 km resolution at 60 degrees north is used. How can I get the right PS for rasnc after which I will apply projectRaster to rasnc to get latlon coordinates

Solution

  • The problem is that rasnc reports that it has a lon/lat coordinate reference system (#coord. ref. : +proj=longlat +datum=WGS84); which looks incorrect. If you know what it should be, (you say UTM, but what is the zone?), assign it:

      crs(rasnc) <- '+proj=utm +zone=??? +datum=WGS84'

    and then things should work. The question is, how did this go wrong?

    I get:

    library(raster)
r <- brick("2013_000_CaPA_continental_v2.3-analyse_final_10km_6h_APCP.nc")
#Loading required namespace: ncdf4
#r
#class       : RasterBrick 
#dimensions  : 824, 935, 770440, 1460  (nrow, ncol, ncell, nlayers)
#resolution  : 10000, 10000  (x, y)
#extent      : -5000, 9345000, -5000, 8235000  (xmin, xmax, ymin, ymax)
#coord. ref. : NA

    which is correct as this file does not provide coordinate reference system information.

    First update your version of raster (you are using ncdf instead of ncdf4, so we can see it is old); and do not leave out code (if you set the coord. ref. system yourself)