Reputation: 2146
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 rasterrize
but could not get the original dimensions of a
. I would like ras
to 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.
Upvotes: 3
Views: 9724
Reputation: 1131
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)
Upvotes: 4