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

[Robert J. Hijmans]

Agus,

If this is from a standard camera, perhaps the problem is that there
is no associated georeference, and different programs may come up with
a different solution to put them somewhere on the map.
R (via gdal) does this: 0, 2640, 0, 1760 (xmin, xmax, ymin, ymax),
with resolution 1 (cell). What does QGIS do?
Perhaps you need to first save it to a geotiff, or add a (fake
coordinates) "world file" so that the treatment becomes consistent.

What is extent(calibf2) ?

Robert

On Thu, Mar 24, 2011 at 12:56 PM, Agustin Lobo <alobolistas at gmail.com> wrote:
> 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
>>>>
>>>
>>
>