A script to mosaic LiDAR-derived DTM and DSMs from the Environment Agency using GRASS.

mosaic_ea_lidar.py

#!/usr/bin/env python2.7
# -*- coding: utf-8 -*-
"""
Mosaic files within a zip archive using GRASS through ARSF DEM Scripts

Designed to mosaic Environemnt Agency LiDAR-derived DTM and DSM files
downloaded in ASCII format from https://environment.data.gov.uk/ds/survey

Usage:

   mosaic_ea_lidar.py -o LIDAR-DSM-2M-SX45.tif LIDAR-DSM-2M-SX45.zip 

Requires arsf_dem_scripts (https://github.com/pmlrsg/arsf_dem_scripts),
GRASS and GDAL Python bindings. See arsf_dem_scripts install instructions.

zip_to_gdal_path function adapted from TuiView (http://tuiview.org/)

Author: Dan Clewley
Creation Date: 02/09/2015

"""
###########################################################
# This file has been created by ARSF Data Analysis Node and
# is licensed under the GPL v3 Licence. A copy of this
# licence is available to download from
# http://www.gnu.org/licenses/gpl-3.0.en.html
###########################################################

from __future__ import print_function
import argparse
import sys
import zipfile
from osgeo import gdal

# Import DEM library
try:
   from arsf_dem import dem_common
   from arsf_dem import dem_utilities
except ImportError as err:
   print("Could not import ARSF DEM library. "
         "You need to install from https://github.com/pmlrsg/arsf_dem_scripts "
         " and make sure it is available on $PYTHONPATH",
         file=sys.stderr)
   print(err, file=sys.stderr)
   sys.exit(1)

EA_LIDAR_NODATA = -9999

def zip_to_gdal_path(filepath):
   """
   Takes in a zip filepath and if the zip contains files
   ascii files, prepend "/viszip" to the path
   so that they can be opened using GDAL without extraction.
   """
   zip_file_list = []

   if zipfile.is_zipfile(filepath):
      try:
         zip_file = zipfile.ZipFile(filepath)
         zip_file_contents = ["/vsizip/{0}/{1}".format(filepath,
                                                zip_info_object.filename) \
                              for zip_info_object in zip_file.filelist \
                              if zip_info_object.filename.endswith(".asc")]
         zip_file_list.extend(zip_file_contents)
         zip_file.close()
      except zipfile.BadZipfile:
         pass

   return zip_file_list

def set_nodata_val(gdal_filename, nodata_val=EA_LIDAR_NODATA):
   """
   Sets No-data value for GDAL dataset
   """
   ds = gdal.Open(gdal_filename, gdal.GA_Update)
   ds.GetRasterBand(1).SetNoDataValue(nodata_val)
   ds = None

if __name__ == "__main__":
   parser = argparse.ArgumentParser(description=\
                                    "A script to mosaic LiDAR-derived DTM and "
                                    "DSM files produced by the EA. "
                                    "Created by ARSF-DAN at "
                                    "Plymouth Marine Laboratory. "
                                    "Latest version available from "
                                    "https://github.com/pmlrsg/arsf_tools/.")
   parser.add_argument("inputfiles", nargs="*",
                       type=str, help="Input zip file(s)")
   parser.add_argument("-o","--outmosaic",
                       type=str, required=True, help="Output mosaic")
   args=parser.parse_args()

   file_list = []

   for zip_file in args.inputfiles:
      file_list.extend(zip_to_gdal_path(zip_file))

   if len(file_list) == 0:
      print("No '.asc' files found within zip file(s) provided as input",
            file=sys.stderr)
      sys.exit(1)

   print("\nCreating mosaic...")
   # Patch files using GRASS
   dem_utilities.patch_files(file_list,
                             args.outmosaic,
                             nodata=EA_LIDAR_NODATA,
                             out_raster_type=dem_common.GDAL_OUTFILE_DATATYPE,
                             projection="UKBNG")
   # Set no-data values (so files display properly)
   set_nodata_val(args.outmosaic)
   print("Finished")