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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
