Converting raster stack extent from meters to lat/lon coordinates

Question

I have a raster stack in the following format which I converted to NetCDF using this method. This works but the latitude and longitude variables are in 'meters' which is the extent of the raster file, and I need them to be in decimal degrees. The R Grid file is here and the format is as follows:

class       : RasterBrick 
dimensions  : 205, 170, 34850, 12  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 160.5, 330.5, 145.5, 350.5  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : \TavgM_1981.grd 
names       :   Jan.1981,   Feb.1981,   Mar.1981,   Apr.1981,   May.1981,   Jun.1981,   Jul.1981,   Aug.1981,   Sep.1981,   Oct.1981,   Nov.1981,   Dec.1981 
min values  :  1.9912137,  0.8120775,  4.0446638,  4.4135274,  6.5349769,  8.6149150,  9.9991916, 11.8400562,  9.6407796,  3.5005649,  4.1205872, -0.6106244 
max values  :   9.850221,   9.121176,  10.238524,   9.858942,  11.669445,  14.260988,  15.722292,  17.235090,  15.708690,  11.598482,  11.552235,   8.981533 
time        : Jan 1981, Feb 1981, Mar 1981, Apr 1981, May 1981, Jun 1981, Jul 1981, Aug 1981, Sep 1981, Oct 1981, Nov 1981, Dec 1981

Ideally I would like to convert this to a NetCDF file which has the following format to be read by a model:

Dimensions:    (lat: 408, lon: 881, nb2: 2, time: 660)
Coordinates:
  * lon        (lon) float64 -44.81 -44.69 -44.56 -44.44 -44.31 -44.19 ...
  * lat        (lat) float64 21.81 21.94 22.06 22.19 22.31 22.44 22.56 22.69 ...
  * time       (time) datetime64[ns] 1951-01-16T12:00:00 1951-02-16T12:00:00 ...
Dimensions without coordinates: nb2
Data variables:
    time_bnds  (time, nb2) datetime64[ns] 1951-01-16T12:00:00 ...
    tas        (time, lat, lon) float64 nan nan nan nan nan nan nan nan nan ...

I tried to convert the Raster Brick extent to lat/lon using this method:

sputm <- SpatialPoints(test@extent@xmin, proj4string=CRS("+proj=utm +zone=29N +datum=WGS84"))

But this gives this error:

Error in (function (classes, fdef, mtable)  : unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric"’

EDIT: Also tried

crs(test) <- '+proj=merc +datum=WGS84'
x <- projectRaster(test, crs='+proj=longlat +datum=WGS84') 

writeRaster(x, "rstack2.nc", overwrite=TRUE, format="CDF",     varname="Temperature", varunit="degC", 
        longname="Temperature -- raster stack to netCDF", xname="Longitude",   yname="Latitude", zname="Time (Month)")

But when I examine this in panoply the lat and lon coordinates are all 0.0.

I'm not sure where to go from here.

The dimensions of your original RasterBrick is strange. Is this really a 1*1m resolution ? In your edit, you force the original CRS to be mercator. How do you know this is the case ? Don't you have any metadata with your original file to know the real CRS ? I cannot use the downloaded file as it lacks the ".grd" file. You can also save it as ".tif" so that we can test it. — Sébastien Rochette, Apr 27 '18 at 06:34
By the way, your "EDIT" is the way to do it. You only need to know the real original CRS — Sébastien Rochette, Apr 27 '18 at 06:35
Thank you - I have contacted the data source to get the original CRS - the only metadata I have with the file is what is printed above (dimensions and extent) - hopefully I will be able to get more. The data should be a 2.5km horizontal grid. I have added the [.grd file here](http://www.filedropper.com/tavgm1981_2) - but i don't think it's possible to get the CRS from this? Maybe it is? — Pad, Apr 27 '18 at 08:26

Sébastien Rochette · Accepted Answer · 2018-04-27T14:52:39.850

3

As fas as I recognized Ireland, I looked at Ireland projections. I think this is the TM75 / Irish Grid. I modified the resolution of your raster so that it is 2500m.

library(raster)

r <- stack("~/Bureau/TavgM_1981")
xmax(r) <- xmin(r) + 2500 * ncol(r)
ymax(r) <- ymin(r) + 2500 * nrow(r)

crs(r) <- "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +ellps=mod_airy +towgs84=482.5,-130.6,564.6,-1.042,-0.214,-0.631,8.15 +units=m +no_defs"

# Convert to degrees
x_wgs84 <- projectRaster(r, crs='+proj=longlat +datum=WGS84') 

writeRaster(x, "~/Bureau/TavgM_wgs84.tif", overwrite = TRUE)

You will see below the result of this resolution + projection modification compared to the target.
We are not so far from reality but there is something bad in your original raster. If this is (unfortunately) usual not to have the crs when somebody gives you a layer, this is not usual to miss the correct resolution.
As the raster you provide is saved as ".gri/.grd" file, this means that it has been produced with R. You should try to get the original R code that produced the file, instead of the final modified layer.

edited Apr 27 '18 at 14:52

answered Apr 27 '18 at 13:48

Sébastien Rochette

6,536
2
22
43

Thank you! I will work on getting the original R code but that may not be possible - I have followed the steps you took but when I examine the file the latitude and longitude are still in meters - I need them to be in decimal degrees - is this possible to do this? I have tried `x <- projectRaster(test, crs='+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +ellps=mod_airy +towgs84=482.5,-130.6,564.6,-1.042,-0.214,-0.631,8.15 +units=m +no_defs') ` but this does not seem to convert correctly – Pad Apr 27 '18 at 14:20
1

I edited the code. Use `projectRaster(r, crs='+proj=longlat +datum=WGS84')` to convert to degrees. However, you now that coordinates will be false. – Sébastien Rochette Apr 27 '18 at 14:54
Thank you so much for all of your help! – Pad Apr 30 '18 at 08:31

Converting raster stack extent from meters to lat/lon coordinates

1 Answers1