2

My question is about using R to transform raster data from a rotated pole projection to a regular latlon grid - and then plotting it.

I am working with data from a numerical weather model, which is specified in a rotated latlon grid - I do this entirely in R. The data is provided in grib format, which is easily read with "raster" package. As example data, I have a grib file in a format identical to the weather data but with ground elevation values for my area of interest (Denmark and southern Sweden) - (example data available here).

Loading and plotting the data looks like:

library("raster")

geopotential <- raster(dk_geopotential_file)

plot(geopotential,
     xlab="Longtitude (rotated)",ylab="Lattitude (rotated)")

Elevation data in the rotated latlon grid.

This looks good - shapes of Denmark and Sweden is clearly seen. Asking for the crs information (as a proj4 string) of the rotated grid gives:

> crs(geopotential)
CRS arguments:
 +proj=longlat +a=6367470 +b=6367470 +no_defs

Now I want to transform this into a regular WGS84 lat-lon format of the type:

"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

I would normally do this with the projectRaster function:

geopotential_wgs84 <- projectRaster(geopotential,
                                    crs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")

But nothing happens! The data remains in the rotated latlon grid. I have tried following these two threads (here and here), and played around with "ob_tran" arguments (Oblique Transformation).

Choosing CRS / PROJ4 string for EURO-Cordex rotated pole projection? seems to do this for a shapefile with the sf package, but I can't make it work with my raster data.

In my case, the coordinates of the geographic South Pole (not the North Pole as in the two linked threads) in the rotated grid is: (lon;lat) = (3;-40).

0

It looks like you will need to include the transformation in the reprojection. Did you happen to see this post? They finally got it in the end using R:

eCRS for rotated pole projection works with rgdal, but not with sf #651

  • Yes, I have been following that post as well. But as I understand, they use vector data and not raster data there? The closest I have gotten to something reasonable is by setting crs(geopotential) <- "+proj=ob_tran +o_proj=longlat +o_lon_p=-177 +o_lat_p=40 +lon_0=180 +to_meter=0.01745329" but the final output from the projectRaster-call still isn't near the regular lonlat coordinates for Denmark – jonaswp Nov 29 '18 at 8:59
0

There is some useful information on Rotation of Coordinates Based On CORDEX Domains and you can contact to AgriMetSoft to get help. Cheers

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.