[R] projectRaster function no values
Adrienne Wootten
amwootte at ncsu.edu
Wed May 18 16:00:49 CEST 2016
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.
> t=1
> 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
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
--
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
On Wed, May 18, 2016 at 9:05 AM, Ben Tupper <btupper at bigelow.org> wrote:
> 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
>
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list