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

[Matteo Mattiuzzi]

Hi,
Robert is right, obeyRes is not required as the function res() does exactly what obeyRes was tought for...
I wrote that function because I had to create rasters with a given resolution to fit as good as possible into a given extent and I was trying to do that with raster() that has no "res" argument + obeyRes. So obeyRes was calculation the available ncol/nrow argument in raster()!
But rasterExtent + res() does the job perfectly!
Thanks Matteo




>>> "Robert J. Hijmans"  01.03.13 19.05 Uhr >>>
There is some confusion here:
To get a google map in a sinusoidial projection you can do

library(dismo)
# extent
e <- extent(c(-1.87710, -1.772096, 55.61399, 55.682171))
# google map
g <- gmap(e, type='hybrid')
#projection
prj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181
+units=m +no_defs"
# project
s1 <- projectRaster(g, crs=prj)
# done



# for more control
# create empty raster
r <- raster(ext=e,crs="+proj=longlat +ellps=WGS84 +datum=WGS84")
# project extent
r <- projectExtent(r,crs=prj)
# set resolution to desired value
# your value of 1000 seemed too large here.
# I do not think you should obeyRes here
# (I am not sure what it is useful for)
res(r) <- 100
# project
s2 <- projectRaster(g, r)


Hope this clarifies.
Robert

On Fri, Mar 1, 2013 at 5:56 AM, Ross Ahmed  wrote:

> 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 at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

    [[alternative HTML version deleted]]

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