Creating a shapefile for footprint of raster's valid data areas with GDAL

Question

I would like to show the borders of the areas of a given ASCII DTM file for which there is data. This problem has been discussed here previously but the recommended solution, the Image Boundary plugin for QGIS, is now considered deprecated. Furthermore, I need to be able to script it, so using the GDAL tools from the command line is probably a better plan for me.

I also have an additional restriction: I need to be able to serve this file on GeoServer, and ultimately display it using Leaflet, so a shapefile or GeoJSON could be good choices for the file format of the final output.

What command should I use?

Answer 1

I have a GRASS/R workflow below that involves GDAL. Its a lot more involved and messy than I recall, but it works so whatever. Blog post here with extra commentary, plain instructions for a Windows environment with OSGeo4W64, R 3.3.1 and GRASS 7.0.5 below.

1: Get all your target rasters into one folder. Make sure that:

• files are named in a logical sequence that does not start with a number.
• files are all in the same projection
• files should be spatially aligned; i.e. the pixels should be all on the same overarching grid system.

2: Go to System > Advanced System Settings > Environment Variables and set the following (adjust as necessary if you installed elsewhere or have an older copy of GRASS):

• GDAL_DATA: C:\OSGeo4W64\share\gdal
• GISBASE: C:\OSGeo4W64\apps\grass\grass-7.0.5
• GISRC: C:\Users\YOUR LOGON HERE\AppData\Roaming\GRASS7\rc
• GRASS_PYTHON: C:\OSGeo4W64\bin\pythonw.exe
• PROJ_LIB: C:\OSGeo4W64\share\proj
• PYTHONPATH: C:\OSGeo4W64\apps\grass\grass-7.0.5\etc\python

• PATH: add the following paths, in the same order as below

  • C:\OSGeo4W64\bin
  • C:\OSGeo4W64\apps\grass\grass-7.0.5\bin
  • C:\OSGeo4W64\apps\grass\grass-7.0.5\lib
  • C:\OSGeo4W64\apps\qgis
  • C:\OSGeo4W64\apps\Python27

• Remove GRASS_ADDON_BASE.

3: Set up GRASS location

Open GRASS and establish a new location with the same CRS as your target rasters. Don’t worry about region settings for now, and there’s no need to create any extra mapsets beyond ‘PERMANENT’. Open the PERMANENT mapset in your new location, and then just exit GRASS. This will update the file that GISRC points to.

R code:

Load libraries and list your target files

library(sp)  
library(raster)  
library(rgdal)  
library(XML)  
library(rgrass7)  

setwd('C:/footprint_test')  
rasters <- list.files(getwd(), pattern = '\\.tif$', full.names = T)  

Import data

grassnames_list <- list()  
for (item in rasters) {  
  grass_name <- substr(basename(item), 1, nchar(basename(item)) - 4)  
  execGRASS('r.in.gdal',   
            flags = 'overwrite',  
            parameters = list(input  = item,  
                              output = grass_name))  
  grassnames_list[[item]] <- grass_name  
}  

extend region to cover imported data

execGRASS('g.region', parameters= list(raster=paste(grassnames_list, collapse=",")))

write an instruction file for r.mapcalc

exp_list <- list()
fpname_list <- list()
for (item in rasters) {
  grassname <- substr(basename(item), 1, nchar(basename(item)) -4)
  fp_name   <- paste0(grassname, '_foot')
  exp       <- paste0(fp_name, ' = if(isnull(', grassname, '), null() , 1)')

  exp_list[[item]]    <- exp
  fpname_list[[item]] <- fp_name
}  

write(c(paste(exp_list, collapse = '\n'), 'end'),
      file = paste0(getwd(), '/rmc_instructions.txt'),
      ncolumns = 1)

use the instruction file to make data presence/absence footprint rasters

rmc   <- 'C:/OSGeo4W64/apps/grass/grass-7.0.5/bin/r.mapcalc.exe'
calcs <- 'C:/footprint_test/rmc_instructions.txt'
system2(rmc, paste0('file=', calcs, ' --overwrite'))

vectorise the footprints

vname_list <- list()
for (fp in fpname_list) {
  v_name <- paste0(fp, '_v')
  execGRASS('r.to.vect', 
            flags = 'overwrite',
            parameters = list(input  = fp,
                              output = v_name,
                              type   = 'area'))
 vname_list[[fp]] <- v_name  
}

write the footprints to file (JSON is also supported)

outname_list <- list()
for (vn in vname_list) {
  out_name <- paste0(getwd(), '/', vn, '.shp')
  execGRASS('v.out.ogr', 
            flags = 'overwrite', 
            parameters = list(input  = vn,
                              output = out_name,
                              format = 'ESRI_Shapefile'))
  outname_list[[vn]] <- out_name
}

If you want multipart polygons

dissname_list <- list()
for (vn in vname_list) {
  diss_name <- paste0(vn, '_sp')
  execGRASS('v.dissolve', 
            flags = 'overwrite',
            parameters = list(input  = vn,
                              column = 'value',
                              output = diss_name))
  dissname_list[[vn]] <- diss_name
}

outname_list2 <- list()
for (dn in dissname_list) {
  out_name2 <- paste0(getwd(), '/', dn, '.shp')
  execGRASS('v.out.ogr', 
            flags = c('m', 'overwrite'), 
            parameters = list(input  = dn,
                              output = out_name2,
                              format = 'ESRI_Shapefile'))
  outname_list2[[dn]] <- out_name2
}

Let me know how you go with that.