Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A script to mosaic LiDAR-derived DTM and DSMs from the Environment Agency using GRASS.
#!/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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.