Skip to content

projectRaster function no values

3 messages · Adrienne Wootten, Ben Tupper

#
All,

Greetings! Any help with this problem is appreciated!

I'm working to get a netcdf file that has a Lambert Conformal Conic
projection into geographic, but also a smaller area.  Here's the issue I'm
having - essentially it looks like projectRaster is working, but the
resulting raster has no values.

The data itself is massive so I can't include that, but here's what's going
on.
time slice of my netcdf
class       : RasterLayer
band        : 1  (of  4  bands)
dimensions  : 229, 234, 53586  (nrow, ncol, ncell)
resolution  : 15, 15  (x, y)
extent      : -1747.5, 1762.5, -1710, 1725  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lon_0=-77 +lat_0=38 +lat_1=30 +lat_2=60
+ellps=WGS84
data source : /projdata/dcerp/DSdata/regcmdata/SE/test.nc
names       : Avg.Max.Aneom.Temperature
z-value     : 1960-01-01
zvar        : TAMAX
Avg.Max.Aneom.Temperature
Min.                   -23.347107
1st Qu.                 -4.706635
Median                   4.347733
3rd Qu.                 16.109032
Max.                    24.556152
NA's                     0.000000
raster with new grid I want
25.94592, 26.08106, 26.21619, 26.35133, 26.48646, 26.62160, 26.75673,
26.89187, 27.02700, 27.16214, 27.29727, 27.43241, 27.56754, 27.70268,
27.83781, 27.97295, 28.10808, 28.24322, 28.37835, 28.51349, 28.64862,
28.78376, 28.91889, 29.05403, 29.18916, 29.32430, 29.45943, 29.59457,
29.72970, 29.86484, 29.99997, 30.13511, 30.27024, 30.40538, 30.54051,
30.67565, 30.81078, 30.94592, 31.08105, 31.21619, 31.35132, 31.48646,
31.62159, 31.75673, 31.89186, 32.02700, 32.16213, 32.29727, 32.43240,
32.56754, 32.70267, 32.83781, 32.97294, 33.10808, 33.24321, 33.37835,
33.51348, 33.64862, 33.78375, 33.91889, 34.05402, 34.18916, 34.32429,
34.45943, 34.59456, 34.72970, 34.86483, 34.99997, 35.13510, 35.27024,
35.40537, 35.54051, 35.67564, 35.81078, 35.94591, 36.08105, 36.21618,
36.35132, 36.48645, 36.62159, 36.75672, 36.89186, 37.02699, 37.16213,
37.29726, 37.43240, 37.56753, 37.70267, 37.83780, 37.97294, 38.10807,
38.24321, 38.37834, 38.51348, 38.64861, 38.78375, 38.91888, 39.05402,
39.18915, 39.32429, 39.45942, 39.59456, 39.72969, 39.86483, 39.99996,
40.13510, 40.27023) # the new regular grid in geographic I'd like to work
with
-102.29720, -102.16207, -102.02693, -101.89180, -101.75666, -101.62153,
-101.48639, -101.35126, -101.21612, -101.08099, -100.94585, -100.81072,
-100.67558, -100.54045, -100.40531, -100.27018, -100.13504,  -99.99991,
 -99.86477, -99.72964, -99.59450, -99.45937, -99.32423, -99.18910,
-99.05396, -98.91883, -98.78369, -98.64856, -98.51342, -98.37829,
-98.24315, -98.10802, -97.97288, -97.83775, -97.70261, -97.56748,
-97.43234, -97.29721, -97.16207, -97.02694, -96.89180, -96.75667,
-96.62153, -96.48640, -96.35126, -96.21613, -96.08099, -95.94586,
-95.81072, -95.67559, -95.54045, -95.40532, -95.27018, -95.13505,
-94.99991, -94.86478, -94.72964, -94.59451, -94.45937, -94.32424,
-94.18910, -94.05397, -93.91883, -93.78370, -93.64856, -93.51343,
-93.37829, -93.24316, -93.10802, -92.97289, -92.83775, -92.70262,
-92.56748, -92.43235, -92.29721, -92.16208, -92.02694, -91.89181,
-91.75667, -91.62154, -91.48640, -91.35127, -91.21613, -91.08100,
-90.94586, -90.81073, -90.67559, -90.54046, -90.40532, -90.27019,
-90.13505, -89.99992, -89.86478, -89.72965, -89.59451, -89.45938,
-89.32424, -89.18911, -89.05397, -88.91884, -88.78370, -88.64857,
-88.51343, -88.37830, -88.24316, -88.10803, -87.97289, -87.83776,
-87.70262, -87.56749, -87.43235, -87.29722, -87.16208, -87.02695,
-86.89181, -86.75668, -86.62154, -86.48641, -86.35127, -86.21614,
-86.08100, -85.94587, -85.81073, -85.67560, -85.54046, -85.40533,
-85.27019, -85.13506, -84.99992, -84.86479, -84.72965, -84.59452,
-84.45938, -84.32425, -84.18911, -84.05398, -83.91884, -83.78371,
-83.64857, -83.51344, -83.37830, -83.24317, -83.10803, -82.97290,
-82.83776, -82.70263, -82.56749, -82.43236, -82.29722, -82.16209,
-82.02695, -81.89182, -81.75668, -81.62155, -81.48641, -81.35128,
-81.21614, -81.08101, -80.94587, -80.81074, -80.67560, -80.54047,
-80.40533, -80.27020, -80.13506, -79.99993, -79.86479, -79.72966,
-79.59452, -79.45939, -79.32425, -79.18912, -79.05398, -78.91885,
-78.78371, -78.64858, -78.51344, -78.37831, -78.24317, -78.10804,
-77.97290, -77.83777, -77.70263, -77.56750, -77.43236, -77.29723,
-77.16209, -77.02696, -76.89182, -76.75669, -76.62155, -76.48642,
-76.35128, -76.21615, -76.08101, -75.94588, -75.81074, -75.67561,
-75.54047, -75.40534, -75.27020, -75.13507, -74.99993, -74.86480,
-74.72966, -74.59453, -74.45939, -74.32426, -74.18912, -74.05399)
dummy raster put to the right extent
project raster itself
class       : RasterLayer
dimensions  : 113, 215, 24295  (nrow, ncol, ncell)
resolution  : 0.1345065, 0.1339391  (x, y)
extent      : -102.9729, -74.05399, 25.13511, 40.27023  (xmin, xmax, ymin,
ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Avg.Max.Aneom.Temperature
values      : NA, NA  (min, max)


I'm quite perplexed with this one, I feel like I'm doing everything right
so I'm not sure what's failing.  The R version is R 3.2.3 in a Linux/Unix
environment.

Many thanks for your help!

Adrienne
#
Hi Adrienne,

You'll always get great help for these kinds of questions if you subscribe and post to the R-SIG-Geo mailing list https://stat.ethz.ch/mailman/listinfo/r-sig-geo  I highly recommend that you join!

I have encountered this before, and usually it is because I have mistakenly assumed that source and destination data are roughly coincident.  I'm not sure if that is true in your case. I have tried to replicate your steps.  I transformed the coordinates of your source and destination rasters into SpatialPoints objects, and then I reprojected your source coordinates to the projection of your destination coordinates.  Unless I have messed up a step, you can see in the plot generated that there is a significant difference in the extent of the source and destination rasters.  Perhaps your destination coordinates are amiss?

Cheers,
Ben


### Start
library(sp)
library(raster)

nc <- 234
nr <- 229
m <- matrix(runif(nc*nr), ncol = nc, nrow = nr)
r1 <- raster(m, xmn = -1747.5, xmx = 1762.5, ymn = -1710, ymx = 1725,
    crs = CRS('+proj=lcc +lon_0=-77 +lat_0=38 +lat_1=30 +lat_2=60'))

newlon <- c(-102.97288, -74.05399)
newlat <- c(25.13511, 40.27023)
newnr <- 113
newnc <- 215
template <- raster(nrows = newnr, ncol = newnc,
    xmn = newlon[1], xmx = newlon[2],
    ymn = newlat[1], ymx = newlat[2],
    crs = CRS("+proj=longlat +datum=WGS84")) 

r2 <- setValues(template, runif(ncell(template)))    

xy_r1 <- SpatialPoints(coordinates(r1),
    proj4string = CRS('+proj=lcc +lon_0=-77 +lat_0=38 +lat_1=30 +lat_2=60'))
xy_r1_tr <- spTransform(xy_r1, CRS("+proj=longlat +datum=WGS84"))


xy_r2 <- SpatialPoints(coordinates(r2),
    proj4string = CRS("+proj=longlat +datum=WGS84"))

plot(xy_r2, pch = '.', axes = TRUE)
points(xy_r1_tr, pch = 1, col = 'orange')

### END
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
#
Ben,

Thanks so for much for the heads up on R-sig-geo (totally forgot about that
one myself).  Thank you also for the clue, it gave me a poke in the right
direction to figure it out.  I'll share with the group just in case someone
else runs across this.  It was actually something off with the resolution
and extent in the source data.
time slice of my netcdf
class       : RasterLayer
band        : 1  (of  4  bands)
dimensions  : 229, 234, 53586  (nrow, ncol, ncell)
resolution  : 15, 15  (x, y)
extent      : -1747.5, 1762.5, -1710, 1725  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=lcc +lon_0=-77 +lat_0=38 +lat_1=30 +lat_2=60
+ellps=WGS84
data source : /projdata/dcerp/DSdata/regcmdata/SE/test.nc
names       : Avg.Max.Aneom.Temperature
z-value     : 1960-01-01
zvar        : TAMAX

The CRS in the file was giving the resolution as 15, but in reality this is
a 15km resolution dataset and should cover an area from 20N to 50N and
about 100W to 50W.  R was interpreting the resolution as 15m, which caused
it to have the small extent, instead of what it should have been, which is
covering Eastern North America.   I had to mess with the resolution and
extent a bit, but once I did it worked beautifully with projectRaster.

Thanks again!

Adrienne