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 <rossah...@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]] _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
