Plotting satellite image in R via sp package -> grid

Dear,

I would like to plot an grid (subset of a satellite image in Sinusoidal
projection) as an spatialobject in R via the sp package (R version 2.5,
winXP).
 
The following works fine:
sinusgrid <- read.table( see for data sample below , header=TRUE,
sep=",", 	na.strings="NA", dec=".", strip.white=TRUE)
sinusgrid <- as.data.frame(sinusgrid)
coordinates(sinusgrid)=~x+y
proj4string(sinusgrid) = CRS("+proj=sinu +lon_0=0 +x_0=0 +y_0=0
+a=6371000 +b=6371000 +units=m")  # Sinusoidal projection
	
l2 = list("SpatialPolygonsRescale", layout.north.arrow(), offset =
c(13416000,-3939000), scale = 500)
l3 = list("SpatialPolygonsRescale", layout.scale.bar(), offset =
c(13416000,-3940000), scale = 1000, fill=c("transparent","black"))
        spplot(sinusgrid, "ndvi", do.log = TRUE,
key.space=list(x=0.2,y=0.8,corner=c(0,1)), scales=list(draw=T), cuts =
4,
     sp.layout=list(l2,l3) )

But when I try ;
gridded(sinusgrid) <- TRUE

I obtain the following error message:
suggested tolerance minimum: 4.31585196603024e-06Error in
points2grid(points, tolerance) : dimension 1 : coordinate intervals are
not constant

before I can use image(sinusgrid) in order to create a spatial map to
which I could overlay vectorlayers and points.

The data is extracted from a satellite image which was projected in
sinusoidal projection (ProjParams =
(6371007.181000,0,0,0,0,0,0,0,21600,0,1,0,0).

How could this satellite data be projected in R such that other shape
files and points can be overlaid on the map?

Regards and thx,
Jan


### Data set: pixels of an satellite image.
x,y,index
13415234.67,-3935966.551,0.83
13415234.67,-3936198.255,0.83
13415234.67,-3936429.96,0.82
13415234.67,-3936661.665,0.87
13415234.67,-3936893.369,0.87
13415234.67,-3937125.074,0.87
13415234.67,-3937356.779,0.88
13415234.67,-3937588.483,0.86
13415234.67,-3937820.188,0.88
13415234.67,-3938051.893,0.77
13415234.67,-3938283.597,0.77
13415234.67,-3938515.302,0.92
13415234.67,-3938747.006,0.89
13415234.67,-3938978.711,0.89
13415234.67,-3939210.416,0.81
13415234.67,-3939442.12,0.76
13415234.67,-3939673.825,0.76
13415234.67,-3939905.53,0.85
13415234.67,-3940137.234,0.85
13415466.38,-3935966.551,0.83
13415466.38,-3936198.255,0.82
13415466.38,-3936429.96,0.83
13415466.38,-3936661.665,0.83
13415466.38,-3936893.369,0.87
13415466.38,-3937125.074,0.89
13415466.38,-3937356.779,0.80
13415466.38,-3937588.483,0.86
13415466.38,-3937820.188,0.88
13415466.38,-3938051.893,0.87

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