Hi Ross,
consider to do it with a GIS, R is not so practical for such data interaction.
Here is what I understand you are trying to do:
 
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]  

Matteo
 

>>> Ross Ahmed <rossahmed@googlemail.com> 01.03.2013 14:56 >>>
Hi all

I have this raster:

myExtent <- extent(c(-1.87710, -1.772096, 55.61399, 55.682171))
r <- raster(ext=myExtent,crs="+proj=longlat +ellps=WGS84 +datum=WGS84")

Ive projected the raster to an equal area projection (Sinusoidal):

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")

Using the obeyRes function kindly given to me by Matteo Mattiuzzi
(https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20130224/3247ffd4/atta
chment.pl), I made each cell in the raster 1km^2:
  
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")

I've downloaded a google map with long/lat projection:

g <- gmap(myExtent, type='hybrid', interpolate=T, lonlat=T)

I'm now trying to project the google map using Sinusoidal (ie, same as
myRaster1km). I've two functions:

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"))
gg <- projectExtent(g,crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181
+b=6371007.181 +units=m +no_defs")

How can I plot the google map using Sinusoidal? Ultimately, my aim is 1.
plot gg 2. plot the myRaster1km over gg 3. draw polygons in each cell of the
myRaster1km and have R recognise which cell I'm drawing in.

Once again, many thanks for help in advance and apologies for any ignorance.

Ross





[[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

	[[alternative HTML version deleted]]

