20 Feb

Merge shapefiles list using OGR (ogr2ogr)

Tired of using ArcGis to merge a huge number of shapefiles?

Here is a simple way to use ogr2ogr and merge shapefiles, so easy!

 
list_files<-Sys.glob('/path_to_shapefiles/*.shp')
 
#create destination file (copy of the first file to merge)
merged_file<-'/path_to_my_results/merged_shapefile.shp'
system(paste(
  'ogr2ogr',
  merged_file,
  list_files[1]))
 
# update the merged_fille by appending all the shapefiles of the list list_files
for (i in seq(2,length(list_files))){
system(paste('ogr2ogr',
             '-update',
             '-append',
             merged_file,
             list_files[i]
             ))
}
04 Dec

Random sampling

Good Evening, fellow GIS analysts!

Here is a quick function to randomly draw a set of points over your area of interest, say for an accuracy assessment.

First the function:

 
randSamp <- function(n,ext,filename=''){
  if(class(ext)=='Extent'){
    x <- runif(n, ext@xmin, ext@xmax)
    y <- runif(n, ext@ymin, ext@ymax)
    s <- SpatialPoints(list(x,y))
  } else if (class(ext) %in% c('Raster','RasterStack','RasterBrick')){
    x <- runif(2*n, extent(ext)@xmin, extent(ext)@xmax)
    y <- runif(2*n, extent(ext)@ymin, extent(ext)@ymax)
    s <- SpatialPoints(list(x,y))
    proj4string(s) <- crs(ext)
  } else if (class(ext) %in% c('SpatialPolygonsDataFrame','SpatialPolygons')){
    x <- c()
    y <- c()
    id <- c()
    while(length(x)<n){
      xl <- runif(2*n, extent(ext)@xmin, extent(ext)@xmax)
      yl <- runif(2*n, extent(ext)@ymin, extent(ext)@ymax)
      sl <- SpatialPoints(list(xl,yl))
      proj4string(sl) <- proj4string(ext)
      point_in_pol <- over( sl , ext , fn = NULL,returnList = FALSE)
      idl <- which(!apply(point_in_pol, 1, function(x) all(is.na(x))))
      x <- cbind(x,xl[idl])
      y <- cbind(y,yl[idl])
    }
    id<-sample(c(1:length(x)),n)
    s <- SpatialPoints(list(x[id],y[id]))
    proj4string(s) <- proj4string(ext)
  } else {
    print('ext not of class "Extent" nor "SpatialPolygonsDataFrame", "SpatialPolygons", "Raster", "RasterStack","RasterBrick". I cannot work with this')
    break
  }
  if(filename!=''){
    s_points_df <- SpatialPointsDataFrame(s, data=data.frame(ID=1:n))
    shapefile(s_poly_df,filename,overwrite=TRUE)
  }
  return(s)
}

And a example on how to use it:

library(sp)
library(rgdal)
library(raster)
set.seed(1)
dat <- matrix(stats::rnorm(2000), ncol = 2)
ch <- chull(dat)
coords <- dat[c(ch, ch[1]), ]  # closed polygon
ext <- SpatialPolygons(list(Polygons(list(Polygon(coords)), ID=1)))
proj4string(ext) = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
spsamp <- randSamp(10,ext)
 
# Now let's visualize this piece of art
plot(ext)
plot(spsamp,add=T,col='red')

sampling

31 Oct

Add ID field to a shapefile using R

Hi GI(rl)S.
Sometimes useful.

library(rgdal)
 
shp<-readOGR("/homes/blackguru/World_Domination.shp", layer="World_Domination")
 
shp@data$ID<-c(1:length(shp@data[,1]))
 
writeOGR(shp, "/homes/blackguru/World_Domination_ID.shp", layer="World_Domination", driver="ESRI Shapefile")
14 Oct

Plot shapefiles with continuous variables

In today’s post, let’s have a look at the plotting functions for continuous variables in shapefiles. To do that we need four functions contained in three packages:
readOGR (rgdal): to read the shapefile
fortify (rgdal): to convert the SpatialPolygonDataFrame to a data frame understandable by ggplot
join (plyr): to join the attributes to the data frame
ggplot (ggplot2): to display the data

The boring part is the conversion from polygon to a data frame that ggplot can handle. To facilitate that, the prepShp4plot function has been created. The hardest part becomes the choice of color for the legend!

prepShp4plot <- function(inShp){
    inShp$id <- as.numeric(1:NROW(inShp))
    inShp.fort <- fortify(inShp, region = "id")
    countryDf <- join(inShp.fort,inShp@data)
}
 
 
# Load the required packages
library(plyr)
library(rgdal)
library(ggplot2)
library(rworldmap)
 
# FOR CLARITY in this example, let's import a shapefile from the rworldmap package
inShp <- getMap(resolution="coarse")
 
# Load the shapefile of interest
#shpDir <- '/home/whiteguru/path2shp'
#setwd(shpDir)
#inShp <- readOGR(dsn='.','world_boundary')
 
 
# Use a small function that prepares your polygon shapefile to be displayed
countryDf <- prepShp4plot(inShp)
 
 
ggplot(data = countryDf, aes(x = long, y = lat, fill = as.numeric(POP_EST), group = group)) +
  geom_polygon(colour = "black") +
  coord_equal() +
  coord_map()+
  theme(panel.background = element_rect(fill='#EBF4FA', colour='black'))+
  scale_fill_gradientn('Human population\n',
                       limits=c(min(na.omit(countryDf$POP_EST)),max(na.omit(countryDf$POP_EST))),
                       colours=c("red","yellow","darkgreen"),
                       breaks=c(min(na.omit(countryDf$POP_EST)),
                       mean(na.omit(countryDf$POP_EST)),
                       max(na.omit(countryDf$POP_EST))),
                       labels = c('Very Small\n','Small\n','Large\n'))


Rplot01