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

[Agustin Lobo]

But then this is not as in gdal. Gdalinfo states that the bottom left is 0,1760
while R states 0,0
Should not raster do as gdal?

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1760.0)
Upper Right ( 2640.0,    0.0)
Lower Right ( 2640.0, 1760.0)
Center      ( 1320.0,  880.0)
Band 1 Block=2640x1 Type=UInt16, ColorInterp=Red
Band 2 Block=2640x1 Type=UInt16, ColorInterp=Green
Band 3 Block=2640x1 Type=UInt16, ColorInterp=Blue

Agus

2011/3/24 Robert J. Hijmans <r.hijmans at gmail.com>:
> 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
>>>>>
>>>>
>>>
>>
>