21 Jan

Dowload OpenStreetMap (OSM) data and read it in QGIS

Happy New GIS,

Be sure to get the last version of GDAL (Ubuntu installation)

Go to Geofabrik.
Download the compressed .osm.pbf files corresponding to our area (Continent, Country).

Using ogr to convert the pbf file into a SQlite database which can be used directly in QGIS with Add SpatiaLite Layer.

ogr2ogr -f "SQLite" -dsco SPATIALITE=YES -spat xmin ymin xmax ymax GuruBlard_Region.db OSM_continent.pbf

where xmin ymin xmax ymax:
spatial query extents, in the SRS of the source layer(s) (or the one specified with -spat_srs). Only features whose geometry intersects the extents will be selected. The geometries will not be clipped unless -clipsrc is specified

For further information, visit this post of Anita Graser.

25 Aug

How to convert hdf files in tif?

Hi crazyGISfans,

Hdf format may be sometimes difficult to manage. We prefer working with geotif.
How do you do convert hdf in tif in R with gdal_translate?

library(gdalUtils)
library(tools)
 
setwd("/guru-blard/")
 
files<-list.files(, pattern=".hdf$")
band<-1
 
for (file in files){
  gdal_translate(src_dataset=get_subdatasets(file)[band], dst_dataset=paste0(file_path_sans_ext(file),"_B", band,'.tif'))
}

That’s it.

28 Jan

Lazy Raster Processing With Gdal Vrts

Hi there, today an introduction to VRT, a magic tool, by Matthew T. Perry. Post from here.

Lazy Raster Processing With Gdal Vrts.
No, not lazy as in REST :-) … Lazy as in “Lazy evaluation”:

    In computer programming, lazy evaluation is the technique of delaying a computation until the result is required.

Take an example raster processing workflow to go from a bunch of tiled, latlong, GeoTiff digital elevation models to a single shaded relief GeoTiff in projected space:

    1) Merge the tiles together
    2) Reproject the merged DEM (using bilinear or cubic interpolation)
    3) Generate the hillshade from the merged DEM

Simple enough to do with GDAL tools on the command line. Here’s the typical, process-as-you-go implementation:

gdal_merge.py -of GTiff -o srtm_merged.tif srtm_12_*.tif 
gdalwarp -t_srs epsg:3310 -r bilinear -of GTiff srtm_merged.tif srtm_merged_3310.tif 
gdaldem hillshade srtm_merged_3310.tif srtm_merged_3310_shade.tif -of GTiff

Alternately, we can simulate lazy evaluation by using GDAL Virtual Rasters (VRT) to perform the intermediate steps, only outputting the GeoTiff as the final step.

gdalbuildvrt srtm_merged.vrt srtm_12_0*.tif
gdalwarp -t_srs epsg:3310 -r bilinear -of VRT srtm_merged.vrt srtm_merged_3310.vrt 
gdaldem hillshade srtm_merged_3310.vrt srtm_merged_3310_shade2.tif -of GTiff

So what’s the advantage to doing it the VRT way? They both produce exactly the same output raster. Lets compare:


benchmark

The Lazy VRT method delays all the computationally-intensive processing until it is actually required. The intermediate files, instead of containing the raw raster output of the actual computation, are XML files which contain the instructions to get the desired output. This allows GDAL to do all the processing in one step (the final step #3). The total processing time is not significantly different between the two methods but in terms of the productivity of the GIS analyst, the VRT method is superior. Imagine working with datasets 1000x this size with many more steps – having to type the command, wait 2 hours, type the next, etc. would be a waste of human resources versus assembling the instructions into vrts then hitting the final processing step when you leave the office for a long weekend.

Additionaly, the VRT method produces only small intermediate xml files instead of having a potentially huge data management nightmare of shuffling around GB (or TB) of intermediate outputs! Plus those xml files serve as an excellent piece of metadata which describe the exact processing steps which you can refer to later or adapt to different datasets.

So next time you have a multi-step raster workflow, use the GDAL VRTs to your full advantage – you’ll save yourself time and disk space by being lazy.

03 Dec

Use gdal_calc to calculate index (NDVI, NDWI, …) , add stats with gdalinfo and overview with gdaladdo

His GISguys!

Here is a simple example of how to use the very fast multi-thred gdal_calc to get NDWI with stats, histogram and overviews.

file='my_input_file.tif'
outfile='my_output_file.tif'
 
  system(
    paste('gdal_calc.py',
          '-A',file,
          '--A_band=2',
          '-B',file,
          '--B_band=4',
          paste(' --outfile=', outpufile,sep=''),
          paste(' --calc=','\"','(A.astype(float)-B)/(A.astype(float)+B)','\"',sep=''),
          paste('--type=','\'','Float32','\'',sep=''),
          paste('--co=','\"','COMPRESS=LZW','\"',sep='')))
 
  system(paste('gdalinfo -hist -stats',outpufile))
 
  system(
    paste('gdaladdo',
          outpufile,
          '2 4 8 16')

Here is another version with byte file as output (values ranging from 0 to 255, faster and smaller files):

file='my_input_file.tif'
outfile='my_output_file.tif'
 
  system(
    paste('gdal_calc.py',
          '-A',file,
          '--A_band=2',
          '-B',file,
          '--B_band=4',
          paste(' --outfile=', outpufile,sep=''),
          paste(' --calc=','\"','((((A.astype(float)-B)/(A.astype(float)+B))+1)*128)','\"',sep=''),
          paste('--type=','\'','Byte','\'',sep=''),
          paste('--co=','\"','COMPRESS=LZW','\"',sep='')))
 
  system(paste('gdalinfo -hist -stats',outpufile))
 
  system(
    paste('gdaladdo',
          outpufile,
          '2 4 8 16')
02 Nov

Efficient Zonal Statistics using R and gdal

Holà crazy GIS lovers,

This function allows to compute zonal statistics of a raster using a ‘zonal shapefile’ and to add an attribute to this shapefile with the result of the statistics for each zone.

NEW VERSION

myZonal <- function (x, z, stat, digits = 0, na.rm = TRUE, 
                     ...) { 
  # source: https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017475.html
  library(data.table)
  fun <- match.fun(stat) 
  vals <- getValues(x) 
  zones <- round(getValues(z), digits = digits) 
  rDT <- data.table(vals, z=zones) 
  setkey(rDT, z) 
  rDT[, lapply(.SD, fun, na.rm = TRUE), by=z] 
} 
 
ZonalPipe<- function (zone.in, raster.in, shp.out=NULL, stat){
  require(raster)
  require(rgdal)
  require(plyr)
 
  # Load raster
  r <- stack(raster.in)
  # Load zone shapefile
  shp <- readOGR(zone.in)
  # Project 'zone' shapefile into the same coordinate system than the input raster
  shp <- spTransform(shp, crs(r))
 
  # Add ID field to Shapefile
  shp@data$ID<-c(1:length(shp@data[,1]))
 
  # Crop raster to 'zone' shapefile extent
  r <- crop(r, extent(shp))	
  # Rasterize shapefile
  zone <- rasterize(shp, r, field="ID", dataType = "INT1U") # Change dataType if nrow(shp) > 255 to INT2U or INT4U
 
  # Zonal stats
  Zstat<-data.frame(myZonal(r, zone, stat))
  colnames(Zstat)<-c("ID", paste0(names(r), "_", c(1:(length(Zstat)-1)), "_",stat))
 
  # Merge data in the shapefile and write it
  shp@data <- plyr::join(shp@data, Zstat, by="ID")
 
  if (is.null(shp.out)){
    return(shp)
  }else{
    writeOGR(shp, shp.out, layer= sub("^([^.]*).*", "\\1", basename(zone.in)), driver="ESRI Shapefile")
  }
}
 
zone.in <- "/home/zone.shp" # Shapefile with zone (INPUT)
raster.in <- "/home/NDVI.tif" #or list.files("/home/, pattern=".tif$") # Raster from which the stats have to be computed (INPUT)
shp.out <- "/home/zone_with_Zstat.shp" # Shapefile with zone + Zonal Stat (OUTPUT)
 
ZonalPipe(zone.in, raster.in, shp.out, stat="mean")
 
shp <- ZonalPipe(zone.in, raster.in, stat="mean")

OLD VERSION

UPDATE: The function has been implemented in the Rnightlight R package.

library(raster)
 
myZonal <- function (x, z, stat, digits = 0, na.rm = TRUE, 
                      ...) { 
  # source: https://stat.ethz.ch/pipermail/r-sig-geo/2013-February/017475.html
  library(data.table)
  fun <- match.fun(stat) 
  vals <- getValues(x) 
  zones <- round(getValues(z), digits = digits) 
  rDT <- data.table(vals, z=zones) 
  setkey(rDT, z) 
  rDT[, lapply(.SD, fun, na.rm = TRUE), by=z] 
} 
 
ZonalPipe<- function (path.in.shp, path.in.r, path.out.r, path.out.shp, zone.attribute, stat){
 
  # 1/ Rasterize using GDAL
 
  #Initiate parameter
  r<-stack(path.in.r)
 
  ext<-extent(r)
  ext<-paste(ext[1], ext[3], ext[2], ext[4])
 
  res<-paste(res(r)[1], res(r)[2])
 
  #Gdal_rasterize
  command<-'gdal_rasterize'
  command<-paste(command, "--config GDAL_CACHEMAX 2000") #Speed-up with more cache (avice: max 1/3 of your total RAM)
  command<-paste(command, "-a", zone.attribute) #Identifies an attribute field on the features to be used for a burn in value. The value will be burned into all output bands.
  command<-paste(command, "-te", as.character(ext)) #(GDAL >= 1.8.0) set georeferenced extents. The values must be expressed in georeferenced units. If not specified, the extent of the output file will be the extent of the vector layers.
  command<-paste(command, "-tr", res) #(GDAL >= 1.8.0) set target resolution. The values must be expressed in georeferenced units. Both must be positive values.
  command<-paste(command, path.in.shp)
  command<-paste(command, path.out.r)
 
  system(command)
 
  # 2/ Zonal Stat using myZonal function
  zone<-raster(path.out.r)
 
  Zstat<-data.frame(myZonal(r, zone, stat))
  colnames(Zstat)[2:length(Zstat)]<-paste0("B", c(1:(length(Zstat)-1)), "_",stat)
 
  # 3/ Merge data in the shapefile and write it
  shp<-readOGR(path.in.shp, layer= sub("^([^.]*).*", "\\1", basename(path.in.shp)))
 
  shp@data <- data.frame(shp@data, Zstat[match(shp@data[,zone.attribute], Zstat[, "z"]),])
 
  writeOGR(shp, path.out.shp, layer= sub("^([^.]*).*", "\\1", basename(path.in.shp)), driver="ESRI Shapefile")
}
 
path.in.shp<-"/home/zone.shp"
path.in.r<-"/home/NDVI.tif" #or path.in.r<-list.files("/home/, pattern=".tif$")
path.out.r<-"/home/zone.tif"
path.out.shp<-"/home/zone_withZstat.shp"
zone.attribute<-"ID"
 
ZonalPipe(path.in.shp, path.in.r, path.out.r, path.out.shp, zone.attribute, stat="mean")
 
#With
#path.in.shp: Shapefile with zone (INPUT)
#path.in.r: Raster from which the stats have to be computed (INPUT)
#path.out.r: Path of path.in.shp converted in raster (intermediate OUTPUT)
#path.out.shp: Path of path.in.shp with stat value (OUTPUT)
#zone.attribute: Attribute name of path.in.shp corresponding to the zones (ID, Country...)
#stat: function to summary path.in.r values ("mean", "sum"...)
24 Oct

Use VRT to deal with very large dataset clip and mosaic

Hi GIS planet,

I have 4000 separeted geoTIFF images of 8000 by 8000 pixels and I want to mosaic and clip those images based on polygons from a shapefile (polygons are 20000 by 20000 pixels).
Here is the code I used:

 
 
## 1 Build a VRT
 
system(
  paste('gdalbuildvrt',
        '/path_to_images/LargeImage.vrt',
        '/path_to_images/*/*.tif')
)
 
 
## 2 Clip & Merge the files based on polygon delimations from a shapefile
zones <- readOGR('/path_to_images/Zones.shp','Zones')
 
for (i in 1:length(zones@polygons))
{
  # get the polygon upper left and down right of the polygones
  ulx<-zones@polygons[[i]]@Polygons[[1]]@coords[2,1]
  uly<-zones@polygons[[i]]@Polygons[[1]]@coords[2,2]
  lrx<-zones@polygons[[i]]@Polygons[[1]]@coords[4,1]
  lry<-zones@polygons[[i]]@Polygons[[1]]@coords[4,2]
 
  system(
    paste('gdal_merge.py',
          paste('-o /path_to_images/LargeImage_Zone_',toString(i),'.tif',sep=''),
          '-ul_lr' ,toString(ulx), toString(uly),toString(lrx), toString(lry),
          '/path_to_images/LargeImage.vrt'))
  system(
    paste('gdaladdo',
          paste('/path_to_images/LargeImage_Zone',toString(i),'.tif',sep=''),
          '2 4 8 16'))
}