[R] projectRaster function no values
Ben Tupper
btupper at bigelow.org
Wed May 18 15:05:17 CEST 2016
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
> On May 17, 2016, at 1:02 PM, Adrienne Wootten <amwootte at ncsu.edu> wrote:
>
> 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.
>
>
>> testvar2 = raster("SE/test.nc",band=t,varname="TAMAX") # pulling first
> time slice of my netcdf
>> testvar2
>
> 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
>
>> summary(testvar2) # yes is does have values
> 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
>
>> newlocs <- raster(ncols=length(newlon), nrows=length(newlat)) # dummy
> raster with new grid I want
>> projection(newlocs)=CRS("+proj=longlat +datum=WGS84")
>
>> newlat=c(25.13511, 25.27025, 25.40538, 25.54052, 25.67565, 25.81079,
> 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
>
>> newlon=c(-102.97288, -102.83774, -102.70261, -102.56747, -102.43234,
> -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)
>
>> extent(newlocs) = c(min(newlon),max(newlon),min(newlat),max(newlat)) #
> dummy raster put to the right extent
>
>> testproj2 = projectRaster(from=testvar2,to=newlocs,method="bilinear") #
> project raster itself
>
>> testproj2 # literally no values.
> 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
> --
> Adrienne Wootten
> Ph.D Candidate / Graduate Research Assistant
> State Climate Office of North Carolina
> Department of Marine, Earth and Atmospheric Sciences
> North Carolina State University
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org
More information about the R-help
mailing list