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

[Robert J. Hijmans]

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
>>
>