[R-sig-Geo] raster: extract is empty with polygon

[Agustin Lobo]

Robert,

I have "reprojection on the fly" not activated. The raster
is an image from an standard camera, hence no projection.
The shapefile takes the projection from the project, but the
coordinates should be within the limits of the image, as it has been
digitized on top of the image. Thus I thought it was just a matter of
setting the CRS for the raster in R.

> intersectExtent(calibf2, SGRGBF40)
Error in intersectExtent(calibf2, SGRGBF40) : Invalid extent

which is in agreement with the rest of results, but still in
contradiction with what we see in QGIS

Agus

2011/3/24 Robert J. Hijmans <r.hijmans at gmail.com>:
> Hi Agus,
>
> It would be useful to see
>
> extent(calibf2)
>
> Apparently it does not overlap with SGRGBF40  (see
> intersectExtent(calibf2, SGRGBF40)  )
> because the polys did not plot on top of the raster.
>
> Perhaps QGIS does some on-the-fly projection trick or other such that
> the coordinates from your on-screen digitization do match those of the
> raster?
>
> This will not fix that! :
>> projection(SGRGBF40) = projection(calibf2)
>
> You can perhaps compare with a polygon you draw in R on top of SGRGBF40
>
> plot(SGRGBF40,1)
> p <- drawPoly()
> extract(SGRGBF40, p)
>
> And save p and load into QGIS?
>
> Best, Robert
>
> ( I would also update the raster package (but I do not think that
> matters for this problem).
>
>
> On Thu, Mar 24, 2011 at 11:27 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
>> I forgot to add that raster and vector do overlap in qgis (the vector
>> was actually digitized on top of the raster)
>> Agus
>>
>> 2011/3/24 Agustin Lobo <alobolistas at gmail.com>:
>>> Hi!
>>> I'm trying to calculate the mean values for a set of polygons.
>>> This is what I do:
>>>
>>> require(raster)
>>> require(rgdal)
>>> #SGRGBWBPS125F40
>>> #"read" tif file
>>> SGRGBF40 = brick("/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif")
>>>> SGRGBF40
>>> class       : RasterBrick
>>> filename    : /media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL/SGRGB/SGRGBWBPS125F40.tif
>>> nlayers     : 3
>>> nrow        : 1760
>>> ncol        : 2640
>>> ncell       : 4646400
>>> projection  : NA
>>> min value   : 0 0 0
>>> max value   : 65535 65535 65535
>>> extent      : 0, 2640, 0, 1760  (xmin, xmax, ymin, ymax)
>>> resolution  : 1, 1  (x, y)
>>>> projection(SGRGBF40)
>>> [1] "NA"
>>>
>>> #read the shape file
>>> calibf2 = readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
>>>> projection(calibf2)
>>> [1] "+proj=utm +zone=31 +ellps=intl +units=m +no_defs"
>>>
>>>> projection(SGRGBF40) = projection(calibf2)
>>>
>>> But they do not overlap as shown by
>>>> plot(subset(SGRGBF40,1))
>>>> plot(calibf2,add=T)
>>>
>>> and the extracted object is empty
>>> #Calculate mean values for each polygon:
>>>> v <- extract(SGRGBF40, calibf2,fun=mean,nl=3)
>>>> summary(v)
>>>      Length Class  Mode
>>>  [1,] 0      -none- NULL
>>>  [2,] 0      -none- NULL
>>>  [3,] 0      -none- NULL
>>>  [4,] 0      -none- NULL
>>>  [5,] 0      -none- NULL
>>> etc
>>>
>>>> sessionInfo()
>>> R version 2.12.2 (2011-02-25)
>>> Platform: i486-pc-linux-gnu (32-bit)
>>>
>>> locale:
>>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>>>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>>>  [5] LC_MONETARY=C              LC_MESSAGES=en_US.UTF-8
>>>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
>>>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
>>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>>
>>> attached base packages:
>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>
>>> other attached packages:
>>> [1] rgdal_0.6-31 raster_1.7-8 sp_0.9-72
>>>
>>> loaded via a namespace (and not attached):
>>> [1] grid_2.12.2     lattice_0.19-17 tools_2.12.2
>>>
>>
>