[R-sig-Geo] spTransform changes the class

[Agustin Lobo]

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.hijmans at 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 <mdsumner at 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 <mdsumner at 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 <alobolistas at 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 at stat.math.ethz.ch
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>