Hi Matteo
Yes that¹s almost exactly what I want. However, I need to be able to see
the fields in the google map, as can be seen here:
myExtent <- extent(c(-1.87710, -1.772096, 55.61399, 55.682171))
g <- gmap(myExtent, type='hybrid', interpolate=T, lonlat=T, zoom=13)
plot(g)
My plan is to draw a polygon around the boundary of each field, and
calculate area of each field.
If the google map cannot be projected to Sinusoidal, perhaps I could first
draw the polygons on the WSG84 google map, then re-project each polygon to
Sinusoidal?
Thanks
Ross
From: Matteo Mattiuzzi <matteo.mattiu...@boku.ac.at>
Date: Friday, 1 March 2013 14:42
To: Ross Ahmed <rossah...@googlemail.com>, <r-sig-geo@r-project.org>
Subject: Antw: [R-sig-Geo] project google map to Sinusoidal
myExtent <- extent(c(-1.87710, -1.772096, 55.61399, 55.682171))
r <- raster(ext=myExtent,crs="+proj=longlat +ellps=WGS84 +datum=WGS84")
myRaster <- projectExtent(r,crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371007.181 +b=6371007.181 +units=m +no_defs")
res1km <- obeyRes(extent(myRaster),resolution=sqrt(1)*1000)
myRaster1km <-
raster(ext=res1km$extent,ncol=res1km$cols,nrows=res1km$rows,crs="+proj=sinu
+lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
myRaster1km[] <- runif(ncell(myRaster1km),0,10)
g <- gmap(myExtent, type='hybrid', interpolate=T, lonlat=T)
gg <- projectRaster(g,crs= "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs")
pols <- rasterToPolygons(myRaster1km)
plot(gg)
plot(pols,add=T)
a <- drawPoly()
projection(a) <- projection(gg)
cell <- cellFromPolygon(myRaster1km,a)[[1]]
myRaster1km[cell]
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
