Who's onlineThere are currently 0 users and 2 guests online.
User loginBook navigationNavigationLive Traffic MapNew Publications
|
Fig. 10.5: 100 simulations of DEM showing using a cross-section from West to East.![]() library(maptools) library(gstat) library(geoR) library(rgdal) library(lattice) library(RSAGA) # MGI Balkans zone 6: gk_6 <- "+proj=tmerc +lat_0=0 +lon_0=18 +k=0.9999 +x_0=6500000 +y_0=0 +ellps=bessel +units=m +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824" # ------------------------------------------------------------ # Import of data to R: # ------------------------------------------------------------ download.file("http://www.spatial-analyst.net/DATA/elevations.zip", destfile=paste(getwd(), "elevations.zip", sep="/")) for(j in list(".shp", ".shx", ".dbf")){ fname <- zip.file.extract(file=paste("elevations", j, sep=""), zipname="elevations.zip") file.copy(fname, paste("./elevations", j, sep=""), overwrite=TRUE) } unlink("elevations.zip") list.files(getwd(), recursive=T, full=F) # import the map to R: elevations <- readShapePoints("elevations.shp", proj4string=CRS(gk_6)) str(elevations) names(elevations@data) <- "Z" # ------------------------------------------------------------ # Variogram modelling: # ------------------------------------------------------------ range(elevations$Z) # sub-sample - geoR can not deal with large datasets! sel <- runif(length(elevations@data[[1]]))<0.2 Z.geo <- as.geodata(elevations[sel,"Z"]) # histogram: plot(Z.geo, qt.col=grey(runif(4))) # Variogram modelling (target variable): par(mfrow=c(1,2)) # anisotropy: plot(variog4(Z.geo, max.dist=1000, messages=FALSE), lwd=2) # fit variogram using likfit: Z.svar <- variog(Z.geo, max.dist=1000, messages=FALSE) Z.vgm <- variofit(Z.svar, messages=FALSE, ini=c(var(Z.geo$data), 1000), fix.nugget=T, nugget=0) Z.vgm # confidence bands: env.model <- variog.model.env(Z.geo, obj.var=Z.svar, model=Z.vgm) plot(Z.svar, envelope=env.model); lines(Z.vgm, lwd=2); dev.off() # spacing between contours: bin.Z <- (max(elevations$Z)-min(elevations$Z))*(length(elevations$Z))^(-1/3) round(bin.Z, 0) # ------------------------------------------------------------ # Geostatistical simulations: # ------------------------------------------------------------ # prepare an empty grid: demgrid <- spsample(elevations, type="regular", cellsize=c(30,30)) gridded(demgrid) <- TRUE fullgrid(demgrid) <- TRUE demgrid@grid gridcell <- demgrid@grid@cellsize[1] # locs <- pred_grid(c(demgrid@bbox[1,1]+gridcell/2, demgrid@bbox[1,2]-gridcell/2), c(demgrid@bbox[2,1]+gridcell/2, demgrid@bbox[2,2]-gridcell/2), by=gridcell) # conditional simulations: # Z.sims <- krige.conv(Z.geo, locations=locs, krige=krige.control(obj.m=Z.vgm), output=output.control(n.predictive=1)) # this is computationally very intensive; geoR is simply not fit to work with large data! # copy the values fitted in geoR: Z.ovgm <- vgm(psill=Z.vgm$cov.pars[1], model="Mat", range=Z.vgm$cov.pars[2], nugget=Z.vgm$nugget, kappa=1.2) # *Kappa is artificially set at 1.2! to produce smooth surface! Z.ovgm # conditional simulations in gstat: N.sim <- 100 DEM.sim <- krige(Z~1, elevations, demgrid, Z.ovgm, nmax=30, nsim=N.sim) fullgrid(DEM.sim) <- TRUE # plot 4 realizations: spplot(DEM.sim[1:4], col.regions=grey(seq(0,1,0.025))) # Cross-section at y=5,073,012: cross.s <- data.frame(X=seq(demgrid@bbox[1,1]+gridcell/2, demgrid@bbox[1,2]-gridcell/2, gridcell), Y=rep(5073012, demgrid@grid@cells.dim[1])) coordinates(cross.s) <-~X+Y # proj4string(cross.s) <- elevations@proj4string cross.ov <- overlay(DEM.sim, cross.s) plot(cross.ov@coords[,1], cross.ov@data[[1]], type="l", xlab="X", ylab="Z", col="grey") for(i in 2:N.sim-1){ lines(cross.ov@coords[,1], cross.ov@data[[i]], col="grey") } lines(cross.ov@coords[,1], cross.ov@data[[N.sim]], lwd=2)
|
Testimonials"From a period in which geographic information systems, and later geocomputation and geographical information science, have been agenda setters, there seems to be interest in trying things out, in expressing ideas in code, and in encouraging others to apply the coded functions in teaching and applied research settings." Poll |
Recent comments
10 years 4 days ago
10 years 18 weeks ago
10 years 26 weeks ago
10 years 38 weeks ago