[R-sig-Geo] Plotting satellite image in R via sp package -> grid AND overlay shape file

[Roger Bivand]

On Thu, 31 May 2007 Jan.Verbesselt at csiro.au wrote:

> Hi Roger,
> 
>  
> 
> Thanks for the help! It works now via the code below. When you now ask
> for the variable GRID, R says that there is no projection system. How
> could a projection be added?  Thx. I would like to add points,  and
> shape files on the georeferenced - map (with scale bar - north arrow?).
> 
>  
> 
> *This is the message:
> 
>  
> 
> >GRID
> 
> ****
> 
> [549,]  13421722  -3939674
> 
> [550,]  13421722  -3939906
> 
> [551,]  13421722  -3940137
> 
> Coordinate Reference System (CRS) arguments: NA
> 
> ***

Thanks, these are longstanding bugs in sp - the projection strings should 
be passed through and are not being. Until we release a new version, 
please do:

proj4string(GRID) <- CRS(
"+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m")

to add it after creation.

> 
>  
> 
> S <- SpatialPoints(cbind( sinusgrid[,1], sinusgrid[,2]), proj4string =
> CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m"))
> 
> SP <- SpatialPixels(S, tolerance =4.3159264735194e-05, proj4string =
> CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m"))
> 
> GRID <- SpatialPixelsDataFrame(points = SP, tolerance
> =4.3159264735194e-05, data = sinusgrid["index"], proj4string =
> CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m"))
> 
> ***
> 
> With this message -->>
> 
> Warning messages:
> 
> 1: grid has empty column/rows in dimension 1 in: points2grid(points,
> tolerance) 
> 
> 2: grid has empty column/rows in dimension 2 in: points2grid(points,
> tolerance)
> 
> ***
> 
>  
> 
> xyz <- as.image.SpatialGridDataFrame(GRID)
> 
> image(GRID, axes=TRUE) # good plot, with coordinates in the axes but
> without a legend! ;/ 
> 
> contour(xyz, add=TRUE)
> 
>  
> 
> In order to add a part of the shape file on top of the map I tried the
> following: (although the grid is only a small part of the shape file!).
> It would be handy to see where in the shape file the grid is positioned.
> 
>  
> 
> library(maptools)
> 
> # --> overlay shape file on top of map (extend of shape file are larger
> than grid)
> 
> pdir <- "C:\\DATA\\shapefile\\Greenhills\\";
> 
> shape<- readShapePoly(paste(pdir,"greenhills_SINUS.shp",sep=""),
> proj4string = CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371000
> +b=6371000 +units=m")) 
> 
> plot(shape) # plot looks fine)
> 
>  
> 
> l3 = list("SpatialPolygonsRescale", layout.scale.bar(), offset
> =c(data[sy,1],data[sy,2]), scale = 500, fill=c("transparent","black")) #
> add scale bar
> 
> rv = list("sp.polygons",shape) # add shape file
> 
> spplot(GRID, c("ndvi"), sp.layout=list(l3,rv))
> 
>  
> 
> * Could the shape file first be plotted and than the grid on top of it? 

Perhaps, but getting the layers as wanted in spplot is not easy.

> 
> * Can a projection be added to the grid - see GRID? 

Do you mean a grid representing some coordinate reference system plotted 
on top of the graphic? In principle yes, but more details would be useful.

> 
> (which other function can be used to fine-tune the creation of spatial
> maps from a combination of a grid and shape file with legend &
> scalebar?)
> 

Knowledge of the underlying graphics engines is necessary, because the 
details are at that level.

Hope this helps,

Roger

>  
> 
> Thanks a lot for your help,
> 
> Jan
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no