## ----setup, include=FALSE, cache=FALSE---------------------------------------- # set global chunk options opts_knit$set(aliases=c(h = 'fig.height', w = 'fig.width'), out.format='latex', use.highlight=TRUE) opts_chunk$set(fig.path='kgraph/fig-', fig.align='center', fig.show='hold', prompt=TRUE, comment=NA, fig.width=6, fig.height=6, out.width='.75\\linewidth', size='scriptsize', message=FALSE) # options to be read by formatR options(replace.assign=TRUE,width=66) ## ----echo=FALSE, results='hide'----------------------------------------------- # make a clean workspace, so packages will show when they load rm(list=ls()) while ("package:dataframes2xls" %in% search()) detach("package:dataframes2xls") while ("package:rgeos" %in% search()) detach("package:rgeos") while ("package:RSQLite" %in% search()) detach("package:RSQLite") while ("package:DBI" %in% search()) detach("package:DBI") while ("package:RColorBrewer" %in% search()) detach("package:RColorBrewer") while ("package:maptools" %in% search()) detach("package:maptools") while ("package:mapdata" %in% search()) detach("package:mapdata") while ("package:maps" %in% search()) detach("package:maps") while ("package:rgdal" %in% search()) detach("package:rgdal") while ("package:raster" %in% search()) detach("package:raster") while ("package:sp" %in% search()) detach("package:sp") ## ----warning=FALSE------------------------------------------------------------ require(sp) require(raster) hwsd <- raster("./HWSD_RASTER/hwsd.bil") ## ----------------------------------------------------------------------------- ncol(hwsd); nrow(hwsd); res(hwsd); extent(hwsd); projection(hwsd) ## ----warning=FALSE------------------------------------------------------------ require(rgdal) (proj4string(hwsd) <-"+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0") ## ----------------------------------------------------------------------------- hwsd.zhnj <- crop(hwsd, extent(c(117.5, 119.5, 31, 33))) nrow(hwsd.zhnj); ncol(hwsd.zhnj); bbox(hwsd.zhnj) ## ----------------------------------------------------------------------------- unique(hwsd.zhnj) ## ----------------------------------------------------------------------------- plot(hwsd.zhnj, col=bpy.colors(length(unique(hwsd.zhnj)))) ## ----warning=FALSE------------------------------------------------------------ hwsd.zhnj3 <- (hwsd.zhnj%/%100) freq(hwsd.zhnj3) require(RColorBrewer) plot(hwsd.zhnj3, col=brewer.pal(length(unique(hwsd.zhnj3)),"Accent")) ## ----find-utm, tidy=FALSE----------------------------------------------------- print(paste("UTM zone:", utm.zone <- floor((sum(bbox(hwsd.zhnj3)[1, ])/2 + 180)/6) + 1)) proj4string.utm50 <- paste("+proj=utm +zone=", utm.zone, "+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0", sep="") hwsd.zhnj3.utm <- projectRaster(hwsd.zhnj3, crs=proj4string.utm50, method="ngb") unique(hwsd.zhnj3.utm) (cell.dim <- res(hwsd.zhnj3.utm)) paste("Cell N dimension is ", round(((cell.dim[2]/cell.dim[1]) - 1)*100,1), "% larger than cell dimension E", sep="") plot(hwsd.zhnj3.utm, col=brewer.pal(6,"Accent"), asp=1) grid() ## ----------------------------------------------------------------------------- (cell.area <- cell.dim[1]*cell.dim[2]/10^4) (tmp <- cbind(freq(hwsd.zhnj3.utm)[,1],freq(hwsd.zhnj3.utm)[,2]*cell.area/10^2)) ix <- which(is.na(tmp[,1])) sum(tmp[-ix,2]) rm(cell.dim, cell.area, tmp, ix) ## ----------------------------------------------------------------------------- rm(hwsd.zhnj3.utm) ## ----echo=FALSE,eval=FALSE---------------------------------------------------- ## # zoom in to find the right pixel for sure ## # do this interactively while writing the doc ## tmp <- crop(hwsd.zhnj, extent(c(118.75, 119, 32, 32.25))) ## plot(tmp) ## (xy <- click(tmp, n=1, id=TRUE, xy=TRUE, type="p")) ## rm(tmp) ## ----eval=FALSE--------------------------------------------------------------- ## plot(hwsd.zhnj, col=bpy.colors(length(unique(hwsd.zhnj)))) ## xy <- click(hwsd.zhnj, n=1, id=TRUE, xy=TRUE, type="p") ## ----echo=FALSE, results='hide'----------------------------------------------- # provide xy as if clicked on Purple Mountain xy <- c(x=(118 + 50/60 + 30/3600), y=(32 + 4/60 + 20/3600), hwsd=11376) xy <- c(x=(118.8292), y=(32.0875), hwsd=11376) ## ----warning=FALSE------------------------------------------------------------ print(xy) (zjs.id <- xy["hwsd"]) rm(xy) ## ----warning=FALSE, tidy=FALSE------------------------------------------------ require(maps) require(mapdata) str(tmp <- map('worldHires','Bhutan', fill=TRUE, plot=FALSE)) require(maptools) bhutan.boundary <- map2SpatialPolygons(tmp, IDs=tmp$names, proj4string= CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) bbox(bhutan.boundary) class(bhutan.boundary) rm(tmp) ## ----------------------------------------------------------------------------- hwsd.bhutan.box <- crop(hwsd, bbox(bhutan.boundary)) hwsd.bhutan.box.sp <- as(hwsd.bhutan.box, "SpatialPixelsDataFrame") sort(unique(hwsd.bhutan.box.sp$hwsd)) ix <- over(hwsd.bhutan.box.sp, bhutan.boundary) hwsd.bhutan.sp <- hwsd.bhutan.box.sp[!is.na(ix),] (bhutan.id <- sort(unique(hwsd.bhutan.sp$hwsd))) ## ----plot.bhutan, out.width='0.48\\linewidth', tidy=FALSE--------------------- spplot(hwsd.bhutan.box.sp, main="HWSD, Bhutan bounding box", col.regions=topo.colors(64), scales=list(draw = TRUE)) spplot(hwsd.bhutan.sp, main="HWSD, Bhutan", col.regions=topo.colors(64), scales=list(draw = TRUE)) ## ----------------------------------------------------------------------------- hwsd.bhutan <- as(hwsd.bhutan.sp, "RasterLayer") rm(hwsd.bhutan.box, hwsd.bhutan.box.sp, ix) ## ----warning=FALSE------------------------------------------------------------ require(RSQLite) m <- dbDriver("SQLite") con <- dbConnect(m, dbname="HWSD.sqlite") dbListTables(con) ## ----------------------------------------------------------------------------- dbGetQuery(con, "pragma table_info(HWSD_DATA)")$name dbGetQuery(con, "pragma table_info(HWSD_DATA)")$type ## ----------------------------------------------------------------------------- dbGetQuery(con, "select count(*) as grid_total from HWSD_DATA") ## ----tidy=FALSE--------------------------------------------------------------- (display.fields <- c("ID","MU_GLOBAL","ISSOIL","SHARE","SU_CODE90", "SU_SYM90","T_USDA_TEX_CLASS")) tmp <- dbGetQuery(con, paste("select", paste(display.fields, collapse=", "), "from HWSD_DATA limit 10")) ## ----------------------------------------------------------------------------- dim(tmp) print(tmp[,display.fields]) rm(tmp) ## ----------------------------------------------------------------------------- str(dbGetQuery(con, "select * from D_SYMBOL90")) ## ----tidy=FALSE--------------------------------------------------------------- (tuple <- dbGetQuery(con, paste("select * from HWSD_DATA where MU_GLOBAL = ", zjs.id))) tuple$SU_SYM90 ## ----tidy=FALSE--------------------------------------------------------------- dbGetQuery(con, paste("select * from D_SYMBOL90 where symbol='", tuple$SU_SYM90,"'",sep="")) ## ----new.zhnj.table, tidy=FALSE----------------------------------------------- dbWriteTable(con, name="WINDOW_ZHNJ", value=data.frame(smu_id=unique(hwsd.zhnj)), overwrite=TRUE) dbGetQuery(con, "pragma table_info(WINDOW_ZHNJ)") ## ----zhnj.join, tidy=FALSE---------------------------------------------------- records <- dbGetQuery(con, "select T.* from HWSD_DATA as T join WINDOW_ZHNJ as U on T.MU_GLOBAL=U.SMU_ID order by SU_SYM90") dim(records) head(records)[,display.fields] sort(unique(records$SU_SYM90)) ## ----------------------------------------------------------------------------- unique(records$SHARE) ## ----------------------------------------------------------------------------- for (i in names(records)[c(2:5,8:15,17:19,28,45)]) { eval(parse(text=paste("records$",i," <- as.factor(records$",i,")", sep=""))) } ## ----------------------------------------------------------------------------- ix <- which(names(records)=="MU_SOURCE2") all(is.na(records[,ix])) ## ----------------------------------------------------------------------------- df <- records for (i in 1:length(names(records))) { if (all(is.na(records[,i]))) df <- df[-i] } dim(records); dim(df) records <- df rm(df, ix, i) ## ----eval=FALSE--------------------------------------------------------------- ## write.csv(records, file="./HWSD_Nanjing.csv") ## ----echo=FALSE, results='hide'----------------------------------------------- write.csv(records, file="./export/HWSD_Nanjing.csv") ## ----eval=FALSE--------------------------------------------------------------- ## require(dataframes2xls) ## write.xls(records, file="./HWSD_Nanjing.xls") ## ----echo=FALSE, results='hide'----------------------------------------------- require(dataframes2xls) write.xls(records, file="./export/HWSD_Nanjing.xls") ## ----see.map.unit.names, tidy=FALSE------------------------------------------- dbExecute(con, "create table TMP as select * from HWSD_DATA as T join WINDOW_ZHNJ as U on T.MU_GLOBAL=U.SMU_ID order by SU_SYM90") head(window.fao90 <- dbGetQuery(con, "select CODE, VALUE, SYMBOL from D_SYMBOL90 as U join TMP as T on T.SU_CODE90=U.CODE")) dbRemoveTable(con, "TMP") ## ----------------------------------------------------------------------------- hwsd.zhnj.sp <- as(hwsd.zhnj, "SpatialGridDataFrame") str(hwsd.zhnj.sp@data) m <- match(hwsd.zhnj.sp@data$hwsd, records$MU_GLOBAL) str(m) hwsd.zhnj.sp@data <- records[m,] rm(m) ## ----plot.nj.raster, out.width='0.48\\linewidth', tidy=FALSE------------------ spplot(hwsd.zhnj.sp, zcol="SU_SYM90", col.regions=topo.colors(length(levels(hwsd.zhnj.sp$SU_SYM90))), main="FAO 1990 soil type code", scales=list(draw = TRUE)) spplot(hwsd.zhnj.sp, zcol="T_SAND", col.regions=bpy.colors(64), main="Toposil sand proportion, %", scales=list(draw = TRUE)) ## ----plot.nj.poly, warning=FALSE, tidy=FALSE---------------------------------- require(rgeos) system.time(hwsd.zhnj.poly <- rasterToPolygons(hwsd.zhnj, n=4, na.rm=TRUE, dissolve=TRUE)) class(hwsd.zhnj.poly) str(hwsd.zhnj.poly@data) spplot(hwsd.zhnj.poly, col.regions=terrain.colors(64), main="HWSD soil map unit ID", scales=list(draw = TRUE)) ## ----------------------------------------------------------------------------- proj4string.utm50 hwsd.zhnj.poly.utm <- spTransform(hwsd.zhnj.poly, CRS(proj4string.utm50)) ## ----------------------------------------------------------------------------- m <- match(hwsd.zhnj.poly.utm$hwsd, records$MU_GLOBAL) str(m) hwsd.zhnj.poly.utm@data <- records[m,] ## ----plot.nj.attributes, out.width='0.48\\linewidth', tidy=FALSE-------------- spplot(hwsd.zhnj.poly.utm, zcol="SU_SYM90", col.regions=topo.colors(length(levels(hwsd.zhnj.sp$SU_SYM90))), main="HWSD FAO 1990 soil type code", scales=list(draw = TRUE)) spplot(hwsd.zhnj.poly.utm, zcol="T_SAND", col.regions=bpy.colors(64), main="Toposil sand proportion, %", scales=list(draw = TRUE)) ## ----bhutan.join, tidy=FALSE-------------------------------------------------- dbWriteTable(con, name="WINDOW_BHUTAN", value=data.frame(smu_id=unique(bhutan.id)), overwrite=TRUE) records.bhutan <- dbGetQuery(con, "select T.* from HWSD_DATA as T join WINDOW_BHUTAN as U on T.MU_GLOBAL=U.SMU_ID order by SU_SYM90") dim(records.bhutan) unique(records.bhutan$MU_GLOBAL) unique(records.bhutan$SHARE) records.bhutan[,c("ID","MU_GLOBAL","SHARE","SU_SYM90","SU_SYM74","T_SAND")] ## ----tidy=FALSE--------------------------------------------------------------- for (i in names(records.bhutan)[c(2:5,8:15,17:19,28,45)]) { eval(parse(text=paste("records.bhutan$",i, " <- as.factor(records.bhutan$",i,")", sep=""))) } ix <- which(names(records.bhutan)=="MU_SOURCE2") df <- records.bhutan for (i in 1:length(names(records.bhutan))) { if (all(is.na(records.bhutan[,i]))) df <- df[-i] } dim(records.bhutan); dim(df) records.bhutan <- df rm(df, ix, i) ## ----eval=FALSE--------------------------------------------------------------- ## write.csv(records.bhutan, file="./HWSD_Bhutan.csv") ## write.xls(records.bhutan, file="./HWSD_Bhutan.xls") ## ----echo=FALSE, results='hide'----------------------------------------------- write.csv(records.bhutan, file="./export/HWSD_Bhutan.csv") write.xls(records.bhutan, file="./export/HWSD_Bhutan.xls") ## ----------------------------------------------------------------------------- hwsd.bhutan.sp <- as(hwsd.bhutan, "SpatialGridDataFrame") m <- match(hwsd.bhutan.sp@data$hwsd, records.bhutan$MU_GLOBAL) hwsd.bhutan.sp@data <- records.bhutan[m, ] summary(hwsd.bhutan.sp@data$T_SAND) ## ----bhutan.tsand, tidy=FALSE------------------------------------------------- spplot(hwsd.bhutan.sp, zcol="T_SAND", col.regions=bpy.colors(64), main="Toposil sand proportion, %, dominant soil", scales=list(draw = TRUE), at=seq(0, 100, 5)) ## ----tidy=FALSE--------------------------------------------------------------- dbExecute(con, "create table TMP as select T.* from HWSD_DATA as T join WINDOW_BHUTAN as U on T.MU_GLOBAL=U.SMU_ID order by SU_SYM90") avg.sand <- dbGetQuery(con, "select MU_GLOBAL, SU_SYM90, SU_SYM74, AVG(T_SAND) as T_SAND_AVG from TMP group by MU_GLOBAL") dim(avg.sand) print(avg.sand) ## ----remove-tmp-avgs---------------------------------------------------------- dbRemoveTable(con, "TMP") ## ----weighted-avgs------------------------------------------------------------ library(dplyr) library(tidyr) bhutan.avg <- records.bhutan %>% dplyr::select(MU_GLOBAL, SHARE, T_SAND, T_CLAY) %>% gather(variable, value, -MU_GLOBAL, -SHARE) %>% mutate(share_2 = ifelse(is.na(value), yes = 0, no = SHARE)) %>% group_by(MU_GLOBAL, variable) %>% mutate(share_2 = share_2 / sum(share_2)) %>% mutate(value = value * share_2) %>% summarise(value = sum(value, na.rm=TRUE)) %>% ungroup() %>% spread(variable, value) print(bhutan.avg) ## ----------------------------------------------------------------------------- hwsd.bhutan.sp <- as(hwsd.bhutan, "SpatialGridDataFrame") m <- match(hwsd.bhutan.sp@data$hwsd, bhutan.avg$MU_GLOBAL) # summary(m) hwsd.bhutan.sp@data <- bhutan.avg[m, ] summary(hwsd.bhutan.sp@data$T_SAND) summary(hwsd.bhutan.sp@data$T_CLAY) ## ----bhutan.tsand.avg, tidy=FALSE--------------------------------------------- spplot(hwsd.bhutan.sp, zcol="T_SAND", col.regions=bpy.colors(64), main="Toposil sand proportion, %, weighted average", scales=list(draw = TRUE), at=seq(0, 100, 5)) ## ----eval=TRUE---------------------------------------------------------------- dbRemoveTable(con, "WINDOW_ZHNJ") dbRemoveTable(con, "WINDOW_BHUTAN") dbDisconnect(con) ## ----echo=FALSE, results='hide'----------------------------------------------- opts_chunk$set(prompt=FALSE) ## ----extract-window, eval=TRUE, tidy=FALSE------------------------------------ ## R script to extract rectangular windows from the Harmonized World Soil Database ## Author: D G Rossiter ## Version: 09-Aug-2017 ##### initialize rm(list=ls()) ##### function to find a UTM zone long2UTM <- function(long) { return(floor((long + 180)/6) + 1) %% 60 } ##### function to extract and format one rectangular window ## arguments: ## bbox: a `raster'-style extent argument, a vector of xmin, xmax, ymin, ymax ## name: a suffix for the file names ## (image, UTM image, csv, excel files, PDF of map unit codes) ## names start with "HWSD_", in subdirectory "window\" and area name ## the image `hwsd' and the SQLite database must ## be already available in the environment extract.one <- function(bbox, name="window") { print(paste("Area name: ", name, "; bounding box: [",paste(bbox,collapse=", "),"]", sep="")) # extract the window dir.create(paste("./window/",name,sep=""), showWarnings = FALSE, recursive=TRUE) setwd(paste("./window/",name,sep="")) hwsd.win <- crop(hwsd, extent(bbox)) # find the zone for the centre of the box print(paste("Central meridian:", centre <- (bbox[1]+bbox[2])/2)) utm.zone <- long2UTM(centre) print(paste("UTM zone:", utm.zone)) # make a UTM version of the window hwsd.win.utm <- projectRaster(hwsd.win, crs=(paste("+proj=utm +zone=",utm.zone, "+datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0", sep="")), method="ngb") print(paste("Cell dimensions:", paste((cell.dim <- res(hwsd.win.utm)), collapse=", "))) # write the raster images to disk eval(parse(text=paste("writeRaster(hwsd.win, file='./HWSD_", name, "', format='EHdr', overwrite=TRUE)",sep=""))) eval(parse(text=paste("writeRaster(hwsd.win.utm, file='./HWSD_", name, "_utm', format='EHdr', overwrite=TRUE)",sep=""))) # extract attributes for just this window dbWriteTable(con, name="WINDOW_TMP", value=data.frame(smu_id=unique(hwsd.win)), overwrite=TRUE) records <- dbGetQuery(con, "select T.* from HWSD_DATA as T join WINDOW_TMP as U on T.mu_global=u.smu_id order by su_sym90") dbRemoveTable(con, "WINDOW_TMP") # convert to factors as appropriate for (i in names(records)[c(2:5,8:15,17:19,28,45)]) { eval(parse(text=paste("records$",i," <- as.factor(records$",i,")", sep=""))) } # remove all-NA fields fields.to.delete <- NULL for (i in 1:length(names(records))) { if (all(is.na(records[,i]))) { fields.to.delete <- c(fields.to.delete, i) } } if (length(fields.to.delete > 1)) records <- records[,-fields.to.delete] print(paste("Dimensions of attribute table: ", paste(dim(records), collapse=", "), " (records, fields with data)", sep="")) # write attribute table in CSV formats eval(parse(text=paste("write.csv(records, file='./HWSD_", name, ".csv')",sep=""))) # make a spatial polygons dataframe, # add attributes print(system.time(hwsd.win.poly <- rasterToPolygons(hwsd.win, n=4, na.rm=TRUE, dissolve=TRUE))) # transform to UTM for correct geometry hwsd.win.poly.utm <- spTransform(hwsd.win.poly, CRS(proj4string(hwsd.win.utm))) m <- match(hwsd.win.poly.utm$value, records$MU_GLOBAL); hwsd.win.poly.utm@data <- records[m,] # plot the map unit ID print(paste("Number of legend categories in the map:", lvls <- length(levels(hwsd.win.poly.utm$MU_GLOBAL)))) eval(parse(text=paste("pdf(file='./HWSD_", name, "_SMU_CODE.pdf')",sep=""))) setwd("../..") } # end extract.one ##### main program ######################################################### ## read in HWSD raster database, assign CRS require(sp) require(raster) hwsd <- raster("./HWSD_RASTER/hwsd.bil") require(rgdal) proj4string(hwsd) <-"+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" ## establish connection to attribute database require(RSQLite) m <- dbDriver("SQLite") con <- dbConnect(m, dbname="HWSD.sqlite") ## other packages to be used in the function require(rgeos) ## call the function for each window we want to extract ### **** NOTE *** change the bounding box: ### c(Long_WestEdge, Long_EastEdge, Lat_SouthEdge, Lat_NorthEdge) ### and also the name of the tile, according to your area ### this example is for the Southern Tier NY/Norther Tier PA (USA) counties extract.one(c(-77, -75, 41, 43), "Twin_Tiers") ## clean up dbDisconnect(con) ## ----extract-country, eval=TRUE, tidy=FALSE----------------------------------- ## R script to extract a county from the Harmonized World Soil Database ## Author: D G Rossiter ## Version: 09-Aug-2017 ##### initialize ######################################################### rm(list=ls()) ##### functions ######################################################### ## function to extract and format one country ## arguments: ## name: a country name, to extract the appropriate bounding polygon(s) ## this name must match the CIA database, see help(worldHires) ## in the `mapdata' package ## will also be used a suffix for the file names ## (image, csv attributes) ## names start with "HWSD_Country_", ## in subdirectory "country\" and area name ## the image `hwsd' and the SQLite database ## must be already available in the environment extract.one <- function(name="") { print(paste("Country:", name)) dir.create(paste("./country/",name,sep=""), showWarnings = FALSE, recursive=TRUE) setwd(paste("./country/",name,sep="")) tmp <- map('worldHires',name, fill=TRUE, plot=FALSE) boundary <- map2SpatialPolygons(tmp, IDs=tmp$names, proj4string= CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) bbox <- bbox(boundary) print(paste("Bounding box: [",paste(t(bbox),collapse=", "),"]", sep="")) # extract the window hwsd.win <- crop(hwsd, extent(bbox)) # overlay only works for sp objects hwsd.win.sp <- as(hwsd.win, "SpatialGridDataFrame") ix <- over(hwsd.win.sp, boundary) hwsd.win.sp <- hwsd.win.sp[!is.na(ix),] hwsd.win <- as(hwsd.win.sp, "RasterLayer") # convert back to raster # find the zone for the centre of the box print(paste("Central meridian:", centre <- (bbox[1]+bbox[2])/2)) # write unprojected raster window image eval(parse(text=paste("writeRaster(hwsd.win, file='./HWSD_raster_", name, "', format='EHdr', overwrite=TRUE)",sep=""))) # extract attributes for this window dbWriteTable(con, name="WINDOW_TMP", value=data.frame(smu_id=unique(hwsd.win)), overwrite=TRUE) records <- dbGetQuery(con, "select T.* from HWSD_DATA as T join WINDOW_TMP as U on T.mu_global=u.smu_id order by su_sym90") dbRemoveTable(con, "WINDOW_TMP") # convert to factors as appropriate for (i in names(records)[c(2:5,8:15,17:19,28,45)]) { eval(parse(text=paste("records$",i," <- as.factor(records$",i,")", sep=""))) } # include all fields print(paste("Dimensions of attribute table: ", paste(dim(records), collapse=", "), " (records, fields)", sep="")) # write attribute table in CSV format eval(parse(text=paste("write.csv(records, file='./HWSD_attributes_", name, ".csv')",sep=""))) setwd("../..") } # end extract.one ##### main program ######################################################### ## read in HWSD raster database, assign CRS require(sp) require(raster) hwsd <- raster("./HWSD_RASTER/hwsd.bil") require(rgdal) proj4string(hwsd) <-"+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" ## establish connection to attribute database require(RSQLite) m <- dbDriver("SQLite") con <- dbConnect(m, dbname="HWSD.sqlite") ## packages for country boundaries require(maps) require(mapdata) require(maptools) ## call the function for each window we want to extract ## *** Note *** replace this with the official name of the country you want ## this name must match the CIA database, see help(worldHires) ## in the `mapdata' package extract.one('Sri Lanka') ## clean up dbDisconnect(con) ## ----eval=FALSE--------------------------------------------------------------- ## (ll <- dbListResults(con)) ## if (length(ll) > 0) { ## res <- ll[[1]] ## if (!dbHasCompleted(res)) dbClearResult(res) ## }