0

I have this raster (less than 10 mb). It is a historical simulation model (SMHI) created by CMIP5 clipped for the city of Amsterdam.

dropbox <- here is the file

What I want to do is to convert this rotated lat/lon to a WGS projection and then re-write the newly obtained raster using ncdf4 or Raster package in R (preferably ncdf4 since I have to do this with many files)

How to do this exactly?

I was trying something like for a single layer:

library(ncdf4)
library(raster)
library(rgdal)

rs <- brick("SMHI_MPI-M-MPI-ESM-LR_historical_r1i1p1_RCA4_amsterdam_1970_2005.nc", varname = 'pr')

rp <- rasterToPoints(rs[[1]])
rp[,3] <-  rp[,3] * 86400 # the precipitation is measured in kg m's (convert here to daily mm)

# change the proj
r2 <-project(rp[,1:2], paste("+proj=ob_tran +o_proj=longlat", "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"), inv=TRUE, use_ob_tran=TRUE)

# add the projected coordinates to the matrix; convert it to df
rp[,1:2] <-r2
rp <-as.data.frame(rp)

# convert it to spatial
coordinates(rp)<-~x+y

# plot it
spplot(rp, scales = list(draw = T))
  • Have you tried creating a raster from your data frame with the projected coords? rps = rasterFromXYZ(as.data.frame(rp))? I get an error reading the NC raster in the first place so I'm not sure it knows about this rotation you mention. – Spacedman Feb 7 at 13:46
  • Hi @Spacedman, unfortunately it doesn't work. It gives an error saying that ' x cell sizes are not regular'. I heard some suggestion saying something to reinterpolate the raster cells however I am not sure. – Andrei Niță Feb 7 at 14:33
  • Hmm I suspect that's the rotation. Can you show a plot of the raster? Does it show the rotation? Possibly related to my "1: In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS: grid_north_pole_latitude=39.25; grid_north_pole_longitude=-162 2" errors when reading into a brick... – Spacedman Feb 7 at 14:42
  • @Spacedman same problem to me. I am able to plot the raster but it does not show the rotation. Only in panoply is possible. – Andrei Niță Feb 7 at 16:05
  • I found something here that I might use. agrimetsoft.com/Cordex%20Coordinate%20Rotation.aspx – Andrei Niță Feb 7 at 16:07

Your Answer

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

Browse other questions tagged or ask your own question.