Search code examples
rmapscoordinatesgisraster

R rasterFromXYZ: x cell sizes are not regular

Given a rasterLayer a,

library(raster)
library(rasterVis)
a=raster('p190001.grd')
head(a)

Based on the description below (see below) which is also found in this link, ftp://ccrp.tor.ec.gc.ca/pub/EC_data/CANGRD/

I came up with the following CRS

mycrs <- CRS("+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-110 +x_0=1884770 +y_0=5220000 +datum=WGS84 +to_meter=50000")

   proj4string(a) <- mycrs
    projection(a) <- mycrs
a
    class       : RasterLayer 
dimensions  : 95, 125, 11875  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -0.5, 124.5, -0.5, 94.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-110 +x_0=1884770 +y_0=5220000 +datum=WGS84 +to_meter=50000 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : p190001 
values      : 4.59, 335.45  (min, max)

a[a==170141000918782798866653488190622531584.00]=NA # change missing value identifier
levelplot(a)

For the sake of completeness, I have created a .Rdata file ~209 KB which can be found here.

The intention is to collapse the a to a dataframe and perform some analysis on it then do a rasterFROMXY. However, I get the error:

Error in rasterFromXYZ(dt) : x cell sizes are not regular

My code:

library(data.table)
v <- getValues(a) # get vector
dt <- as.data.table(v)

# Add LATITUDE and LONGITUDE columns to dt

dt[,LATITUDE:=grid_pnt_lls$y]
dt[,LONGITUDE:=grid_pnt_lls$x]
dtnames <- c(names(dt)[2:3],names(dt)[1]) 

setcolorder(dt,dtnames)

vcols <- c('LONGITUDE','LATITUDE','v')

setcolorder(dt,vcols)

library(data.table)
setnames(dt, "LONGITUDE", "x")
setnames(dt, "LATITUDE", "y")
setnames(dt, "v", "z")


ras=rasterFromXYZ(dt)#     Error in rasterFromXYZ(dt) : x cell sizes are not regular

As a workaround, I tried to projectRaster to regular latlon (projj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0") ) and even rasterrizebut could not get the original dimensions of a. I would like rasto have same dimensions as a but on latlon grid.

#===========================================================================

The CANGRD grid is in polar stereographic projection with a 50 km spatial resolution. The grid is a 125 (columns) by 95 (rows) matrix, where the SW corner (0,0) is at 40.0451°N latitude and 129.8530°W longitude. The projection is true at 60.0°N and centered on 110.0°W. The file ‘CANGRD_points_LL.txt’ lists the latitudes and longitudes for each grid point.

The general format of the ‘YYYYDD.grd’ file is:
Id – ‘DSAA’ identifies the file as an ASCII grid file
nx ny - nx is the integer number of grid lines along the X axis (columns)
        ny is the integer number of grid lines along the Y axis (rows)
xlo xhi - xlo is the minimum X value of the grid
          xhi is the maximum X value of the grid
ylo yhi - ylo is the minimum Y value of the grid
          yhi is the maximum Y value of the grid
zlo zhi - are the minimum and maximum Z values of the grid.
Solution

  • Instead of using getValues use as.data.frame to convert raster to a dataframe. 

    dt <- data.table(as.data.frame(a, xy = TRUE))
setnames(dt, "p190001", "z")

## Sample Calculation on Values
dt[, z := z/2]

ras <- rasterFromXYZ(dt)
plot(ras)

    enter image description here