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.mattiuzzi@boku.ac.at>
Date:  Friday, 1 March 2013 14:42
To:  Ross Ahmed <rossahmed@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]]

