3

I have a set of rasters that I've had trouble reprojecting to Robinson. It fails without errors, the crs changes but the extent doesn't convert to metres and it doesn't plot properly. Sources say EPSG:32663 is deprecated.

Projecting raster with '+init=EPSG:4326' and the proj4 string together doesn't make a difference

r <- projectRaster(r, crs='+init=EPSG:4326 +proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

But, setting the crs this way, followed by projecting with the proj4 string does

crs(r) <- "+init=EPSG:4326"


r <- projectRaster(r, crs='+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

I need to do this across a list of rasters but am having trouble getting the syntax right for setting the CRS inside a function, example

for (i in 1:length(filenames)){
   crs(filenames[[i]]) <- "+init=EPSG:4326"
   # projectRaster(eval(as.name(filenames[[i]])),crs='+init=EPSG:4326 +proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
}

gives the following error

Error in (function (classes, fdef, mtable)  :  unable to find an inherited method for function ‘crs<-’ for signature ‘"character"’

This is the raster I'm working with

ftp://ftp.star.nesdis.noaa.gov/pub/sod/mecb/crw/data/5km/v3.1/nc/v1.0/monthly/2019/ct5km_baa-max_v3.1_201901.nc

class      : RasterLayer 
dimensions : 3600, 7200, 25920000  (nrow, ncol, ncell)
resolution : 0.05, 0.05  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +init=EPSG:32663 +proj=eqc +lat_ts=0 +lat_0=0 
+lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
+ellps=WGS84 +towgs84=0,0,0 
source     : /Users/traveler/Desktop/Reefs/alert- 
animation/raw/ct5km_baa-max_v3.1_201607.nc 
names      : bleaching.alert.area 
z-value    : 1122811200 
zvar       : bleaching_alert_area 
7
  • 2
    Why can't you read the raster in and use projectRaster without assign? What goes wrong with the code you have? Why have you pointed us to one raster when your problem seems to be with several rasters? Have you got several like this? Like that in what way?
    – Spacedman
    Aug 15, 2019 at 13:29
  • 1
    If this file really was in epsg:32663 then its units, and hence its extent should be in metres not degrees. When you project it to Robinson and it seems to fail its because it thinks its a 360 metre x 180 metre piece of the earth centred on (0,0), not a full global map. So whoever wrote this has lost some important data somewhere. If it is epsg:4326 then try overlaying a known good epsg:4326 vector map and check the alignment.
    – Spacedman
    Aug 15, 2019 at 14:26
  • 1
    Hi, I edited my question. If I try to change the CRS in the function I get "Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘crs<-’ for signature ‘"character"’ and if I try to reproject I get no error, but no result either. Here is the rest of the data (monthly files) ftp.star.nesdis.noaa.gov/pub/sod/mecb/crw/data/5km/v3.1/nc/v1.0/… Aug 15, 2019 at 15:34
  • That error is because you are trying to set the crs on the filename - read the file into a raster object then override the incorrect 32663 CRS with what you think is the correct one (4326, but check that).
    – Spacedman
    Aug 15, 2019 at 17:21
  • filenames is a list of my rasters after being read in, I'm passing them to the crs() function - are there other ways to override it? Aug 15, 2019 at 20:30

1 Answer 1

1

Try something like:

readandproject = function(filenames){
 lapply(filenames, function(filename){
   r = raster(filename)
   crs(r) = "+init=epsg:4326" # fix broken CRS
   r = projectRaster(r, crs="+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
   return(r)
   }
)
}

Then prasters = readandproject(c("file1.nc","file2.nc")) should return a list of two projected rasters from the filenames. plot(prasters[[1]]) to see the first one.

Also, I loaded one of the rasters into QGIS and overrode the CRS to WGS84 (epsg:4326) and compared it with OpenStreetMap's base layer, and the correlation looked pretty good up to the resolution of the raster.

2
  • This is elegant - oddly spTransform gives an "unable to find an inherited method for function" error but using projectRaster here enables it to run. I have a country borders shapefile that I have reprojected and added to the plot and it matches up nicely - is this enough to be confident that it is 4326? (I really appreciate your help, I'd still be stuck right now!) Aug 16, 2019 at 8:07
  • Yeah sorry, I was ad-libbing. projectRaster it should be I'll edit.
    – Spacedman
    Aug 16, 2019 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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