[R-sig-Geo] Antw: project google map to Sinusoidal

[Matteo Mattiuzzi]

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]