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")
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)
[1] "NA"
#read the shape file
calibf2 = readOGR(dsn="/media/TRANSCEND/MONTSENY2008/CALIBRACION2010/DigiPreprocess/TESTCASA/CALSEL",layer="calibf2",stringsAsFactors=F)
[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
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