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

[Agustin Lobo]

Anyway, the problem is solved by including the appropriate wld file:
1.0
0
0
-1.0
0.5
0.5

In this case, both QGIS and R (raster) keep the same coordinates and
the vector layer matches the raster layer.

Agus

2011/3/25 Agustin Lobo <alobolistas at gmail.com>:
> 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
>>>>>>
>>>>>
>>>>
>>>
>>
>