Dear OSGEO users,

to settle this topic  - as Roger Bivand, who responded to my query advised 
- I shortly summarize my problems again and put some pieces of code for a 
final solution to it.
The case was to use the output of the gridded values from the spatial 
population genetics package Geneland for the production of a multiple 
layers plot which show some georeferenced vector layers, a map in lon/lat 
coordinates and a overlaid raster with the output results from Geneland. 
The layers overlaid, in this case, was first a partial map of Europe 
covering the sampling area of the example, a couple of overlaid rivers for 
orientation, a spatial points object displaying the sampling points, and at 
last the assignment probabilities which were output of Geneland.
Coming back to my first problem where I was running off of the registration 
of the Geneland output. This happened in the beginning, because the grid in 
an additional step was transformed to a raster ascii (ArcView format) which 
bounding box was not correct.
The solution to this was - as Roger proposed - to directly use the output 
in a SPatialGrid object.
The following piece of code should give an outline of the final solution.
Thanks for the help.
#+#+#+#+#+#+#+#+#+#+#+##+#+#+#+#
require(maptools)
require(rgdal)


file.exists("c:\\GEOSTATISTIK\\ESRI_Daten\\EpiInfo_Europe\\eur_trans\\eur_trans.shp")
?readShapePoly
eur.SPpdf <- 
readShapePoly("c:\\GEOSTATISTIK\\ESRI_Daten\\europe_transalpina\\europa_auswahl.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
plot(eur.SPpdf)
?readShapeLines
eur.river.SPldf <- 
readShapeLines("c:\\AxelHille\\shapefiles\\DEU_rheinmaindonau.shp",
proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 
+no_defs"))        # WGS84
plot(eur.river.SPldf)

#LONLAT
#?rgdal::spTransform
mat.neu.LL
str(assign.tra3.LL)    #[, c(3,3)]
tra3LL.xy.spdf <- SpatialPointsDataFrame(mat.neu.LL, 
assign.tra3.LL[,c(3,3)], proj4string=CRS("+proj=longlat +ellps=WGS84 
+datum=WGS84 +no_defs"))
summary(tra3LL.xy.spdf)
plot(tra3LL.xy.spdf)

################################################################################
GENELAND output

################################################################################
path.mcmc=MCMC.PATH.TRANS3LL
  geneland_proba_df <- read.table(paste(path.mcmc, 
"proba.pop.membership.txt", sep = ""))
  str(geneland_proba_df)
  coords <- geneland_proba_df[, 1:2]
  head(coords)
  assign.prob <- matrix(geneland_proba_df[, 3], nrow=100, ncol=100, byrow 
= FALSE)
  assign.prob
library(sp)
?sp
#get size of the grid in Geneland from parameter file
    param.postprocess <- as.matrix(read.table(paste(path.mcmc, 
        "postprocess.parameters.txt", sep = "")))
    nxdom <- as.numeric(param.postprocess[1, 3])
    nxdom
    nydom <- as.numeric(param.postprocess[2, 3])

#xygrid <- expand.grid(1:100, 1:100) # in this case the Geneland grid was 
100 times 100
# in general
xygrid <- expand.grid(1:nxdom, 1:nydom)

SPDF <- SpatialPointsDataFrame(SpatialPoints(coords,
proj4string=CRS("+proj=longlat +datum=WGS84")),
data=data.frame(z=c(z)))
summary(SPDF)
str(SPDF, maxlevel=1)
bb <- bbox(SPDF)
bb
grd <- GridTopology(bb[,1], c(diff(bb[1,])/nxdom, diff(bb[2,])/nydom), 
c(nxdom,nydom))
grd
SG <- SpatialGrid(grd, CRS("+proj=longlat +datum=WGS84"))
o <- overlay(SG, SPDF)
oo <- tapply(SPDF$z, o, mean, na.rm=TRUE)
oo
df <- data.frame(z=c(z))
#is.data.frame(df)
df$z[as.integer(names(oo))] <- oo
SGDF <- SpatialGridDataFrame(grd, CRS("+proj=longlat +datum=WGS84"), 
data=df)
str(SGDF)
####################
    
-- 
Dr. rer. nat. Axel Hille
Diplom-Biologe

GRATIS für alle GMX-Mitglieder: Die maxdome Movie-FLAT!


	[[alternative HTML version deleted]]

