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

***

 

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? 

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

(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?)

 

Thanks a lot for your help,

Jan


	[[alternative HTML version deleted]]


