JezK
Edit File: ogr_layer_algebra.py
#!/usr/bin/env python3 # -*- coding: utf-8 -*- # ****************************************************************************** # $Id$ # # Project: GDAL Python Interface # Purpose: Application for executing OGR layer algebra operations # Author: Even Rouault, even dot rouault at spatialys.com # # ****************************************************************************** # Copyright (c) 2013, Even Rouault <even dot rouault at spatialys.com> # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING # FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER # DEALINGS IN THE SOFTWARE. # ****************************************************************************** import os import sys from osgeo import gdal, ogr, osr ############################################################################### def Usage(): print( """ Usage: ogr_layer_algebra.py Union|Intersection|SymDifference|Identity|Update|Clip|Erase -input_ds name [-input_lyr name] -method_ds [-method_lyr name] -output_ds name [-output_lyr name] [-overwrite] [-opt NAME=VALUE]* [-f format_name] [-dsco NAME=VALUE]* [-lco NAME=VALUE]* [-input_fields NONE|ALL|fld1,fl2,...fldN] [-method_fields NONE|ALL|fld1,fl2,...fldN] [-nlt geom_type] [-a_srs srs_def]""" ) return 2 ############################################################################### def EQUAL(a, b): return a.lower() == b.lower() ############################################################################### def CreateLayer( output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt, ): output_lyr = output_ds.CreateLayer(output_lyr_name, srs, geom_type, lco) if output_lyr is None: print('Cannot create layer "%s"' % output_lyr_name) return None if input_fields == "ALL" and method_fields == "ALL": return output_lyr input_prefix = "" method_prefix = "" for val in opt: if val.lower().find("input_prefix=") == 0: input_prefix = val[len("input_prefix=") :] elif val.lower().find("method_prefix=") == 0: method_prefix = val[len("method_prefix=") :] if input_fields == "ALL": layer_defn = input_lyr.GetLayerDefn() for idx in range(layer_defn.GetFieldCount()): fld_defn = layer_defn.GetFieldDefn(idx) fld_defn = ogr.FieldDefn( input_prefix + fld_defn.GetName(), fld_defn.GetType() ) if output_lyr.CreateField(fld_defn) != 0: print( 'Cannot create field "%s" in layer "%s"' % (fld_defn.GetName(), output_lyr.GetName()) ) elif input_fields != "NONE": layer_defn = input_lyr.GetLayerDefn() for fld in input_fields: idx = layer_defn.GetFieldIndex(fld) if idx < 0: print( 'Cannot find field "%s" in layer "%s"' % (fld, layer_defn.GetName()) ) continue fld_defn = layer_defn.GetFieldDefn(idx) fld_defn = ogr.FieldDefn( input_prefix + fld_defn.GetName(), fld_defn.GetType() ) if output_lyr.CreateField(fld_defn) != 0: print( 'Cannot create field "%s" in layer "%s"' % (fld, output_lyr.GetName()) ) if method_fields == "ALL": layer_defn = method_lyr.GetLayerDefn() for idx in range(layer_defn.GetFieldCount()): fld_defn = layer_defn.GetFieldDefn(idx) fld_defn = ogr.FieldDefn( method_prefix + fld_defn.GetName(), fld_defn.GetType() ) if output_lyr.CreateField(fld_defn) != 0: print( 'Cannot create field "%s" in layer "%s"' % (fld_defn.GetName(), output_lyr.GetName()) ) elif method_fields != "NONE": layer_defn = method_lyr.GetLayerDefn() for fld in method_fields: idx = layer_defn.GetFieldIndex(fld) if idx < 0: print( 'Cannot find field "%s" in layer "%s"' % (fld, layer_defn.GetName()) ) continue fld_defn = layer_defn.GetFieldDefn(idx) fld_defn = ogr.FieldDefn( method_prefix + fld_defn.GetName(), fld_defn.GetType() ) if output_lyr.CreateField(fld_defn) != 0: print( 'Cannot create field "%s" in layer "%s"' % (fld, output_lyr.GetName()) ) return output_lyr ############################################################################### def main(argv=sys.argv): version_num = int(gdal.VersionInfo("VERSION_NUM")) if version_num < 1100000: print("ERROR: Python bindings of GDAL 1.10 or later required") return 1 driver_name = "ESRI Shapefile" quiet = False input_ds_name = None input_lyr_name = None method_ds_name = None method_lyr_name = None output_ds_name = None output_lyr_name = None op_str = None dsco = [] lco = [] opt = [] overwrite = False input_fields = "ALL" method_fields = None geom_type = ogr.wkbUnknown srs_name = None srs = None argv = ogr.GeneralCmdLineProcessor(argv) if argv is None: return 1 # Parse command line arguments. i = 1 while i < len(argv): arg = argv[i] if arg == "-f" and i + 1 < len(argv): i = i + 1 driver_name = argv[i] elif arg == "-input_ds" and i + 1 < len(argv): i = i + 1 input_ds_name = argv[i] elif arg == "-input_lyr" and i + 1 < len(argv): i = i + 1 input_lyr_name = argv[i] elif arg == "-method_ds" and i + 1 < len(argv): i = i + 1 method_ds_name = argv[i] elif arg == "-method_lyr" and i + 1 < len(argv): i = i + 1 method_lyr_name = argv[i] elif arg == "-output_ds" and i + 1 < len(argv): i = i + 1 output_ds_name = argv[i] elif arg == "-output_lyr" and i + 1 < len(argv): i = i + 1 output_lyr_name = argv[i] elif arg == "-input_fields" and i + 1 < len(argv): i = i + 1 if EQUAL(argv[i], "NONE"): input_fields = "NONE" elif EQUAL(argv[i], "ALL"): input_fields = "ALL" else: input_fields = argv[i].split(",") elif arg == "-method_fields" and i + 1 < len(argv): i = i + 1 if EQUAL(argv[i], "NONE"): method_fields = "NONE" elif EQUAL(argv[i], "ALL"): method_fields = "ALL" else: method_fields = argv[i].split(",") elif arg == "-dsco" and i + 1 < len(argv): i = i + 1 dsco.append(argv[i]) elif arg == "-lco" and i + 1 < len(argv): i = i + 1 lco.append(argv[i]) elif arg == "-opt" and i + 1 < len(argv): i = i + 1 opt.append(argv[i]) elif arg == "-nlt" and i + 1 < len(argv): i = i + 1 val = argv[i] if EQUAL(val, "NONE"): geom_type = ogr.wkbNone elif EQUAL(val, "GEOMETRY"): geom_type = ogr.wkbUnknown elif EQUAL(val, "POINT"): geom_type = ogr.wkbPoint elif EQUAL(val, "LINESTRING"): geom_type = ogr.wkbLineString elif EQUAL(val, "POLYGON"): geom_type = ogr.wkbPolygon elif EQUAL(val, "GEOMETRYCOLLECTION"): geom_type = ogr.wkbGeometryCollection elif EQUAL(val, "MULTIPOINT"): geom_type = ogr.wkbMultiPoint elif EQUAL(val, "MULTILINESTRING"): geom_type = ogr.wkbMultiLineString elif EQUAL(val, "MULTIPOLYGON"): geom_type = ogr.wkbMultiPolygon elif EQUAL(val, "GEOMETRY25D"): geom_type = ogr.wkbUnknown | ogr.wkb25DBit elif EQUAL(val, "POINT25D"): geom_type = ogr.wkbPoint25D elif EQUAL(val, "LINESTRING25D"): geom_type = ogr.wkbLineString25D elif EQUAL(val, "POLYGON25D"): geom_type = ogr.wkbPolygon25D elif EQUAL(val, "GEOMETRYCOLLECTION25D"): geom_type = ogr.wkbGeometryCollection25D elif EQUAL(val, "MULTIPOINT25D"): geom_type = ogr.wkbMultiPoint25D elif EQUAL(val, "MULTILINESTRING25D"): geom_type = ogr.wkbMultiLineString25D elif EQUAL(val, "MULTIPOLYGON25D"): geom_type = ogr.wkbMultiPolygon25D else: print("-nlt %s: type not recognised." % val) return 1 elif arg == "-a_srs" and i + 1 < len(argv): i = i + 1 srs_name = argv[i] elif EQUAL(arg, "Union"): op_str = "Union" elif EQUAL(arg, "Intersection"): op_str = "Intersection" elif EQUAL(arg, "SymDifference"): op_str = "SymDifference" elif EQUAL(arg, "Identity"): op_str = "Identity" elif EQUAL(arg, "Update"): op_str = "Update" elif EQUAL(arg, "Clip"): op_str = "Clip" elif EQUAL(arg, "Erase"): op_str = "Erase" elif arg == "-overwrite": overwrite = True elif arg == "-q" or arg == "-quiet": quiet = True else: return Usage() i = i + 1 if ( input_ds_name is None or method_ds_name is None or output_ds_name is None or op_str is None ): return Usage() if method_fields is None: if op_str in ("Update", "Clip", "Erase"): method_fields = "NONE" else: method_fields = "ALL" if input_fields == "NONE" and method_fields == "NONE": print( "Warning: -input_fields NONE and -method_fields NONE results in all fields being added" ) # Input layer input_ds = ogr.Open(input_ds_name) if input_ds is None: print("Cannot open input dataset : %s" % input_ds_name) return 1 if input_lyr_name is None: cnt = input_ds.GetLayerCount() if cnt != 1: print( "Input datasource has not a single layer, so you should specify its name with -input_lyr" ) return 1 input_lyr = input_ds.GetLayer(0) else: input_lyr = input_ds.GetLayerByName(input_lyr_name) if input_lyr is None: print('Cannot find input layer "%s"' % input_lyr_name) return 1 # Method layer method_ds = ogr.Open(method_ds_name) if method_ds is None: print("Cannot open method dataset : %s" % method_ds_name) return 1 if method_lyr_name is None: cnt = method_ds.GetLayerCount() if cnt != 1: print( "Method datasource has not a single layer, so you should specify its name with -method_lyr" ) return 1 method_lyr = method_ds.GetLayer(0) else: method_lyr = method_ds.GetLayerByName(method_lyr_name) if method_lyr is None: print('Cannot find method layer "%s"' % method_lyr_name) return 1 # SRS if srs_name is not None: if not EQUAL(srs_name, "NULL") and not EQUAL(srs_name, "NONE"): srs = osr.SpatialReference() if srs.SetFromUserInput(srs_name) != 0: print("Failed to process SRS definition: %s" % srs_name) return 1 else: srs = input_lyr.GetSpatialRef() srs2 = method_lyr.GetSpatialRef() if srs is None and srs2 is not None: print("Warning: input layer has no SRS defined, but method layer has one.") elif srs is not None and srs2 is None: print( "Warning: input layer has a SRS defined, but method layer has not one." ) elif srs is not None and srs2 is not None and srs.IsSame(srs2) != 1: print( "Warning: input and method layers have SRS defined, but they are not identical. No on-the-fly reprojection will be done." ) # Result layer output_ds = ogr.Open(output_ds_name, update=1) if output_ds is None: output_ds = ogr.Open(output_ds_name) if output_ds is not None: print( 'Output datasource "%s" exists, but cannot be opened in update mode' % output_ds_name ) return 1 drv = ogr.GetDriverByName(driver_name) if drv is None: print("Cannot find driver %s" % driver_name) return 1 output_ds = drv.CreateDataSource(output_ds_name, options=dsco) if output_ds is None: print('Cannot create datasource "%s"' % output_ds_name) return 1 # Special case if ( EQUAL(drv.GetName(), "ESRI Shapefile") and output_lyr_name is None and EQUAL(os.path.splitext(output_ds_name)[1], ".SHP") ): output_lyr_name = os.path.splitext(os.path.basename(output_ds_name))[0] if output_lyr_name is None: print("-output_lyr should be specified") return 1 output_lyr = CreateLayer( output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt, ) if output_lyr is None: return 1 else: drv = output_ds.GetDriver() if output_lyr_name is None: cnt = output_ds.GetLayerCount() if cnt != 1: print( "Result datasource has not a single layer, so you should specify its name with -output_lyr" ) return 1 output_lyr = output_ds.GetLayer(0) output_lyr_name = output_lyr.GetName() else: output_lyr = output_ds.GetLayerByName(output_lyr_name) if output_lyr is None: if ( EQUAL(drv.GetName(), "ESRI Shapefile") and EQUAL(os.path.splitext(output_ds_name)[1], ".SHP") and not EQUAL( output_lyr_name, os.path.splitext(os.path.basename(output_ds_name))[0], ) ): print( 'Cannot create layer "%s" in a shapefile called "%s"' % (output_lyr_name, output_ds_name) ) return 1 output_lyr = CreateLayer( output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt, ) if output_lyr is None: return 1 if overwrite: cnt = output_ds.GetLayerCount() iLayer = None # initialize in case there are no loop iterations for iLayer in range(cnt): poLayer = output_ds.GetLayer(iLayer) if poLayer is not None and poLayer.GetName() == output_lyr_name: break if iLayer != cnt: if output_ds.DeleteLayer(iLayer) != 0: print("DeleteLayer() failed when overwrite requested.") return 1 output_lyr = CreateLayer( output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt, ) if output_lyr is None: return 1 op = getattr(input_lyr, op_str) if not quiet: ret = op(method_lyr, output_lyr, options=opt, callback=gdal.TermProgress_nocb) else: ret = op(method_lyr, output_lyr, options=opt) input_ds = None method_ds = None output_ds = None if ret != 0: print("An error occurred during %s operation" % op_str) return 1 return 0 if __name__ == "__main__": sys.exit(main(sys.argv))