Ross, here the area using rgeos...but again I think a GIS System is much more suited for this type of things... Matteo
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, rgb=TRUE, zoom=13) gg <- projectRaster(g, method="ngb",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) plotRGB(gg) plot(pols,add=T) a <- drawPoly() projection(a) <- projection(gg) require(rgeos) gArea(a) # area of a in mapunit [qm] cell <- cellFromPolygon(myRaster1km,a)[[1]] >>> Ross Ahmed 03/01/13 5:09 PM >>> 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 Date: Friday, 1 March 2013 14:42 To: Ross Ahmed , 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] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo