Thank you all for your ideas. Please note though that the goal is to create a raster as a reference object for running resample(). The original data is a high resolution raster in EPSG 3035 (Br) that I want to resample to a lon,lat grid of 0.25 deg resolution (Frankly, I'm astonished that such a simple operation is so complicated in R, I think that it would be very simple in grass)
After having created the grid eugrd025, my first idea was just considering that a grid is not but a set of adjacent rectangular polygons and hence raster::polygonValues() would solve my problem. Thus I converted the grid eugrd025 to SpatialPolygons (eugrd025DF), reprojected to EPSG 3035 (eugrd025EFDC_SPDF) and tried polygonValues() . This did not work (see [1]) and I was advised to try an "all raster approach", which is a good approximation in this case as the resolution of the grid is much larger than the resolution of the original raster Br. Therefore, I needed a raster version of SpPolygons object eugrd025EFDC_SPDF. raster() does not accept SpPolygons (which makes sense) , but SpPixels or SpGrids. Unfortunately, while you can convert from SpPixels to SpPolygons, the reverse, i.e., eugrd025EFDC_SPX <- as(eugrd025EFDC_SPDF, "SpatialPixels") does not work Indefatigable, I decided to use the original grid eugrd025, reproject it via spTransform(), convert to raster using raster() and then resample(), but spTransform() had a low blow ready for me: spTransform() does not work with Spgrids, its help page rgdal/html/spTransform-methods.html does not mention it and it silently converts the output to Sppoints, which ultimately results on a wrong resample() result. Finally, the following mixture of Mike and Robert advices worked (fast) > eugrd025r <- raster(eugrd025) > Br025rEFDC <- projectRaster(Br,eugrd025r) Thanks! Agus [1] http://r-sig-geo.2731867.n2.nabble.com/polygonValues-raster-Very-slow-tp5239776p5239776.html 2010/7/1 Robert J. Hijmans <r.hijm...@gmail.com>: > The general idea is that for changing the projection of a raster you > provide the input RasterLayer and a Raster object with the spatial > parameters (projection, extent, resolution) to which it should be > transformed. For the latter, you can use the RasterLayer that you want > to project the input data to (I call it 'x'). > > prjgrd <- projectRaster(eugrd025, x, method='bilinear') > > Note that it should *not* be necessary to resample after projection > because projecting a raster already implies resampling. > > Robert > > On Thu, Jul 1, 2010 at 8:48 AM, Michael Sumner <mdsum...@gmail.com> wrote: >> Ah, sorry that was completely wrong. I think this is right, but still >> untested: >> >> raster.eugrd025 <- raster(eugrd025) >> >> pr <- projectExtent(raster.eugrd025, projection(eugrd025EFDC_SPDF)) >> >> eugrd025EFDC <- projectRaster(raster.eugrd025, pr) >> >> >> >> On Fri, Jul 2, 2010 at 1:46 AM, Michael Sumner <mdsum...@gmail.com> wrote: >>> spTransform cannot reproject a grid - this will (usually) requires >>> destructive resampling of the data to a completely new grid. There are >>> functions in the raster package to do this: projectRaster. >>> >>> I think this would work, but there may be important details that you >>> will need to investigate: >>> >>> raster.eugrd025 <- raster(eugrd025) >>> >>> eugrd025EFDC <- projectRaster(raster.eugrd025, >>> CRSobj=CRS(projection(eugrd025EFDC_SPDF))) >>> >>> On Fri, Jul 2, 2010 at 1:40 AM, Agustin Lobo <alobolis...@gmail.com> wrote: >>>> Hi! >>>> >>>> I use spTransform to reproject: >>>> >>>>> class(eugrd025) >>>> [1] "SpatialGrid" >>>> attr(,"package") >>>> [1] "sp" >>>> >>>>> projection(eugrd025) >>>> [1] "+proj=longlat +ellps=WGS84" >>>> >>>>> eugrd025EFDC <- spTransform(eugrd025, >>>>> CRSobj=CRS(projection(eugrd025EFDC_SPDF))) >>>> >>>>> projection(eugrd025EFDC) >>>> [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 >>>> +ellps=GRS80 +units=m +no_defs" >>>> >>>> This is good, but: >>>>> class(eugrd025EFDC) >>>> [1] "SpatialPoints" >>>> attr(,"package") >>>> [1] "sp" >>>> >>>> Why is the class changed to SpatialPoints? I need a grid object for >>>>> eugrd025EFDCr <- raster(eugrd025EFDC) >>>>> Br025 <- resample(Br,eugrd025EFDCr) >>>> >>>> don't I? >>>> >>>> I also have the problem that the CRS of eugrd025EFDCr is not conserved: >>>>> projection(Br) >>>> [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 >>>> +ellps=GRS80 +units=m +no_defs" >>>> >>>>> projection(eugrd025EFDC) >>>> [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 >>>> +ellps=GRS80 +units=m +no_defs" >>>> >>>>> projection(eugrd025EFDCr) >>>> [1] "NA" >>>> >>>> which I fix with >>>>> projection(eugrd025EFDCr) <- CRS(projection(eugrd025EFDC)) >>>> >>>> and repeat >>>>> Br025 <- resample(Br,eugrd025EFDCr) >>>> >>>> Anyway, I get a much coarser raster than I want: >>>>> dim(Br025) >>>> [1] 10 10 1 >>>> >>>> because >>>>> dim(eugrd025EFDCr) >>>> [1] 10 10 1 >>>> >>>> While it should be 180x100: >>>> >>>>> summary(eugrd025EFDC) >>>> Object of class SpatialPoints >>>> Coordinates: >>>> min max >>>> s1 2498538 6561151 >>>> s2 1327858 4313263 >>>> Is projected: TRUE >>>> proj4string : >>>> [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 >>>> +units=m +no_defs] >>>> Number of points: 18000 >>>>> dim(Br025) >>>> [1] 10 10 1 >>>> >>>>> summary(eugrd025) >>>> Object of class SpatialGrid >>>> Coordinates: >>>> min max >>>> coords.x1 -10.125 34.875 >>>> coords.x2 34.875 59.875 >>>> Is projected: FALSE >>>> proj4string : [+proj=longlat +ellps=WGS84] >>>> Number of points: 2 >>>> Grid attributes: >>>> cellcentre.offset cellsize cells.dim >>>> 1 -10 0.25 180 >>>> 2 35 0.25 100 >>>> >>>> Is eugrd025EFDCr wrong because class(eugrd025EFDC) has been changed >>>> from SpatialGrid to SpatialPoints? >>>> >>>> Agus >>>> >>>> _______________________________________________ >>>> R-sig-Geo mailing list >>>> R-sig-Geo@stat.math.ethz.ch >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>> >>> >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > _______________________________________________ R-sig-Geo mailing list R-sig-Geo@stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
