1#!/usr/bin/env python3 2# -*- coding: utf-8 -*- 3# ****************************************************************************** 4# $Id: ogr_layer_algebra.py e4fe7cc06270e5f38dfe78e6785a6bcca4e39e29 2021-04-01 21:02:04 +0300 Idan Miara $ 5# 6# Project: GDAL Python Interface 7# Purpose: Application for executing OGR layer algebra operations 8# Author: Even Rouault, even dot rouault at spatialys.com 9# 10# ****************************************************************************** 11# Copyright (c) 2013, Even Rouault <even dot rouault at spatialys.com> 12# 13# Permission is hereby granted, free of charge, to any person obtaining a 14# copy of this software and associated documentation files (the "Software"), 15# to deal in the Software without restriction, including without limitation 16# the rights to use, copy, modify, merge, publish, distribute, sublicense, 17# and/or sell copies of the Software, and to permit persons to whom the 18# Software is furnished to do so, subject to the following conditions: 19# 20# The above copyright notice and this permission notice shall be included 21# in all copies or substantial portions of the Software. 22# 23# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 24# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 25# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 26# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 27# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 28# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 29# DEALINGS IN THE SOFTWARE. 30# ****************************************************************************** 31 32import os 33import sys 34 35from osgeo import gdal 36from osgeo import ogr 37from osgeo import osr 38 39############################################################################### 40 41 42def Usage(): 43 print(""" 44Usage: ogr_layer_algebra.py Union|Intersection|SymDifference|Identity|Update|Clip|Erase 45 -input_ds name [-input_lyr name] 46 -method_ds [-method_lyr name] 47 -output_ds name [-output_lyr name] [-overwrite] 48 [-opt NAME=VALUE]* 49 [-f format_name] [-dsco NAME=VALUE]* [-lco NAME=VALUE]* 50 [-input_fields NONE|ALL|fld1,fl2,...fldN] [-method_fields NONE|ALL|fld1,fl2,...fldN] 51 [-nlt geom_type] [-a_srs srs_def]""") 52 return 1 53 54############################################################################### 55 56 57def EQUAL(a, b): 58 return a.lower() == b.lower() 59 60############################################################################### 61 62 63def CreateLayer(output_ds, output_lyr_name, srs, geom_type, lco, 64 input_lyr, input_fields, 65 method_lyr, method_fields, opt): 66 67 output_lyr = output_ds.CreateLayer(output_lyr_name, srs, geom_type, lco) 68 if output_lyr is None: 69 print('Cannot create layer "%s"' % output_lyr_name) 70 return None 71 72 input_prefix = '' 73 method_prefix = '' 74 for val in opt: 75 if val.lower().find('input_prefix=') == 0: 76 input_prefix = val[len('input_prefix='):] 77 elif val.lower().find('method_prefix=') == 0: 78 method_prefix = val[len('method_prefix='):] 79 80 if input_fields == 'ALL': 81 layer_defn = input_lyr.GetLayerDefn() 82 for idx in range(layer_defn.GetFieldCount()): 83 fld_defn = layer_defn.GetFieldDefn(idx) 84 fld_defn = ogr.FieldDefn(input_prefix + fld_defn.GetName(), fld_defn.GetType()) 85 if output_lyr.CreateField(fld_defn) != 0: 86 print('Cannot create field "%s" in layer "%s"' % (fld_defn.GetName(), output_lyr.GetName())) 87 88 elif input_fields != 'NONE': 89 layer_defn = input_lyr.GetLayerDefn() 90 for fld in input_fields: 91 idx = layer_defn.GetFieldIndex(fld) 92 if idx < 0: 93 print('Cannot find field "%s" in layer "%s"' % (fld, layer_defn.GetName())) 94 continue 95 fld_defn = layer_defn.GetFieldDefn(idx) 96 fld_defn = ogr.FieldDefn(input_prefix + fld_defn.GetName(), fld_defn.GetType()) 97 if output_lyr.CreateField(fld_defn) != 0: 98 print('Cannot create field "%s" in layer "%s"' % (fld, output_lyr.GetName())) 99 100 if method_fields == 'ALL': 101 layer_defn = method_lyr.GetLayerDefn() 102 for idx in range(layer_defn.GetFieldCount()): 103 fld_defn = layer_defn.GetFieldDefn(idx) 104 fld_defn = ogr.FieldDefn(method_prefix + fld_defn.GetName(), fld_defn.GetType()) 105 if output_lyr.CreateField(fld_defn) != 0: 106 print('Cannot create field "%s" in layer "%s"' % (fld_defn.GetName(), output_lyr.GetName())) 107 108 elif method_fields != 'NONE': 109 layer_defn = method_lyr.GetLayerDefn() 110 for fld in method_fields: 111 idx = layer_defn.GetFieldIndex(fld) 112 if idx < 0: 113 print('Cannot find field "%s" in layer "%s"' % (fld, layer_defn.GetName())) 114 continue 115 fld_defn = layer_defn.GetFieldDefn(idx) 116 fld_defn = ogr.FieldDefn(method_prefix + fld_defn.GetName(), fld_defn.GetType()) 117 if output_lyr.CreateField(fld_defn) != 0: 118 print('Cannot create field "%s" in layer "%s"' % (fld, output_lyr.GetName())) 119 120 return output_lyr 121 122############################################################################### 123 124 125def main(argv=None): 126 version_num = int(gdal.VersionInfo('VERSION_NUM')) 127 if version_num < 1100000: 128 print('ERROR: Python bindings of GDAL 1.10 or later required') 129 return 1 130 131 frmt = 'ESRI Shapefile' 132 quiet_flag = 0 133 input_ds_name = None 134 input_lyr_name = None 135 method_ds_name = None 136 method_lyr_name = None 137 output_ds_name = None 138 output_lyr_name = None 139 op_str = None 140 dsco = [] 141 lco = [] 142 opt = [] 143 overwrite = False 144 input_fields = 'ALL' 145 method_fields = None 146 geom_type = ogr.wkbUnknown 147 srs_name = None 148 srs = None 149 150 argv = ogr.GeneralCmdLineProcessor(argv) 151 if argv is None: 152 return 1 153 154 # Parse command line arguments. 155 i = 1 156 while i < len(argv): 157 arg = argv[i] 158 159 if arg == '-f' and i + 1 < len(argv): 160 i = i + 1 161 frmt = argv[i] 162 163 elif arg == '-input_ds' and i + 1 < len(argv): 164 i = i + 1 165 input_ds_name = argv[i] 166 167 elif arg == '-input_lyr' and i + 1 < len(argv): 168 i = i + 1 169 input_lyr_name = argv[i] 170 171 elif arg == '-method_ds' and i + 1 < len(argv): 172 i = i + 1 173 method_ds_name = argv[i] 174 175 elif arg == '-method_lyr' and i + 1 < len(argv): 176 i = i + 1 177 method_lyr_name = argv[i] 178 179 elif arg == '-output_ds' and i + 1 < len(argv): 180 i = i + 1 181 output_ds_name = argv[i] 182 183 elif arg == '-output_lyr' and i + 1 < len(argv): 184 i = i + 1 185 output_lyr_name = argv[i] 186 187 elif arg == '-input_fields' and i + 1 < len(argv): 188 i = i + 1 189 if EQUAL(argv[i], "NONE"): 190 input_fields = "NONE" 191 elif EQUAL(argv[i], "ALL"): 192 input_fields = "ALL" 193 else: 194 input_fields = argv[i].split(',') 195 196 elif arg == '-method_fields' and i + 1 < len(argv): 197 i = i + 1 198 if EQUAL(argv[i], "NONE"): 199 method_fields = "NONE" 200 elif EQUAL(argv[i], "ALL"): 201 method_fields = "ALL" 202 else: 203 method_fields = argv[i].split(',') 204 205 elif arg == '-dsco' and i + 1 < len(argv): 206 i = i + 1 207 dsco.append(argv[i]) 208 209 elif arg == '-lco' and i + 1 < len(argv): 210 i = i + 1 211 lco.append(argv[i]) 212 213 elif arg == '-opt' and i + 1 < len(argv): 214 i = i + 1 215 opt.append(argv[i]) 216 217 elif arg == "-nlt" and i + 1 < len(argv): 218 i = i + 1 219 val = argv[i] 220 221 if EQUAL(val, "NONE"): 222 geom_type = ogr.wkbNone 223 elif EQUAL(val, "GEOMETRY"): 224 geom_type = ogr.wkbUnknown 225 elif EQUAL(val, "POINT"): 226 geom_type = ogr.wkbPoint 227 elif EQUAL(val, "LINESTRING"): 228 geom_type = ogr.wkbLineString 229 elif EQUAL(val, "POLYGON"): 230 geom_type = ogr.wkbPolygon 231 elif EQUAL(val, "GEOMETRYCOLLECTION"): 232 geom_type = ogr.wkbGeometryCollection 233 elif EQUAL(val, "MULTIPOINT"): 234 geom_type = ogr.wkbMultiPoint 235 elif EQUAL(val, "MULTILINESTRING"): 236 geom_type = ogr.wkbMultiLineString 237 elif EQUAL(val, "MULTIPOLYGON"): 238 geom_type = ogr.wkbMultiPolygon 239 elif EQUAL(val, "GEOMETRY25D"): 240 geom_type = ogr.wkbUnknown | ogr.wkb25DBit 241 elif EQUAL(val, "POINT25D"): 242 geom_type = ogr.wkbPoint25D 243 elif EQUAL(val, "LINESTRING25D"): 244 geom_type = ogr.wkbLineString25D 245 elif EQUAL(val, "POLYGON25D"): 246 geom_type = ogr.wkbPolygon25D 247 elif EQUAL(val, "GEOMETRYCOLLECTION25D"): 248 geom_type = ogr.wkbGeometryCollection25D 249 elif EQUAL(val, "MULTIPOINT25D"): 250 geom_type = ogr.wkbMultiPoint25D 251 elif EQUAL(val, "MULTILINESTRING25D"): 252 geom_type = ogr.wkbMultiLineString25D 253 elif EQUAL(val, "MULTIPOLYGON25D"): 254 geom_type = ogr.wkbMultiPolygon25D 255 else: 256 print("-nlt %s: type not recognised." % val) 257 return 1 258 259 elif arg == "-a_srs" and i + 1 < len(argv): 260 i = i + 1 261 srs_name = argv[i] 262 263 elif EQUAL(arg, "Union"): 264 op_str = "Union" 265 266 elif EQUAL(arg, "Intersection"): 267 op_str = "Intersection" 268 269 elif EQUAL(arg, "SymDifference"): 270 op_str = "SymDifference" 271 272 elif EQUAL(arg, "Identity"): 273 op_str = "Identity" 274 275 elif EQUAL(arg, "Update"): 276 op_str = "Update" 277 278 elif EQUAL(arg, "Clip"): 279 op_str = "Clip" 280 281 elif EQUAL(arg, "Erase"): 282 op_str = "Erase" 283 284 elif arg == "-overwrite": 285 overwrite = True 286 287 elif arg == '-q' or arg == '-quiet': 288 quiet_flag = 1 289 290 else: 291 return Usage() 292 293 i = i + 1 294 295 if input_ds_name is None or \ 296 method_ds_name is None or \ 297 output_ds_name is None or \ 298 op_str is None: 299 return Usage() 300 301 if method_fields is None: 302 if op_str in ('Update', 'Clip', 'Erase'): 303 method_fields = 'NONE' 304 else: 305 method_fields = 'ALL' 306 307 if input_fields == 'NONE' and method_fields == 'NONE': 308 print('Warning: -input_fields NONE and -method_fields NONE results in all fields being added') 309 310 # Input layer 311 input_ds = ogr.Open(input_ds_name) 312 if input_ds is None: 313 print('Cannot open input dataset : %s' % input_ds_name) 314 return 1 315 316 if input_lyr_name is None: 317 cnt = input_ds.GetLayerCount() 318 if cnt != 1: 319 print('Input datasource has not a single layer, so you should specify its name with -input_lyr') 320 return 1 321 input_lyr = input_ds.GetLayer(0) 322 else: 323 input_lyr = input_ds.GetLayerByName(input_lyr_name) 324 if input_lyr is None: 325 print('Cannot find input layer "%s"' % input_lyr_name) 326 return 1 327 328 # Method layer 329 method_ds = ogr.Open(method_ds_name) 330 if method_ds is None: 331 print('Cannot open method dataset : %s' % method_ds_name) 332 return 1 333 334 if method_lyr_name is None: 335 cnt = method_ds.GetLayerCount() 336 if cnt != 1: 337 print('Method datasource has not a single layer, so you should specify its name with -method_lyr') 338 return 1 339 method_lyr = method_ds.GetLayer(0) 340 else: 341 method_lyr = method_ds.GetLayerByName(method_lyr_name) 342 if method_lyr is None: 343 print('Cannot find method layer "%s"' % method_lyr_name) 344 return 1 345 346 # SRS 347 if srs_name is not None: 348 if not EQUAL(srs_name, "NULL") and not EQUAL(srs_name, "NONE"): 349 srs = osr.SpatialReference() 350 if srs.SetFromUserInput(srs_name) != 0: 351 print("Failed to process SRS definition: %s" % srs_name) 352 return 1 353 else: 354 srs = input_lyr.GetSpatialRef() 355 srs2 = method_lyr.GetSpatialRef() 356 if srs is None and srs2 is not None: 357 print('Warning: input layer has no SRS defined, but method layer has one.') 358 elif srs is not None and srs2 is None: 359 print('Warning: input layer has a SRS defined, but method layer has not one.') 360 elif srs is not None and srs2 is not None and srs.IsSame(srs2) != 1: 361 print('Warning: input and method layers have SRS defined, but they are not identical. No on-the-fly reprojection will be done.') 362 363 # Result layer 364 output_ds = ogr.Open(output_ds_name, update=1) 365 if output_ds is None: 366 output_ds = ogr.Open(output_ds_name) 367 if output_ds is not None: 368 print('Output datasource "%s" exists, but cannot be opened in update mode' % output_ds_name) 369 return 1 370 371 drv = ogr.GetDriverByName(frmt) 372 if drv is None: 373 print('Cannot find driver %s' % frmt) 374 return 1 375 376 output_ds = drv.CreateDataSource(output_ds_name, options=dsco) 377 if output_ds is None: 378 print('Cannot create datasource "%s"' % output_ds_name) 379 return 1 380 381 # Special case 382 if EQUAL(drv.GetName(), "ESRI Shapefile") and output_lyr_name is None \ 383 and EQUAL(os.path.splitext(output_ds_name)[1], ".SHP"): 384 output_lyr_name = os.path.splitext(os.path.basename(output_ds_name))[0] 385 386 if output_lyr_name is None: 387 print('-output_lyr should be specified') 388 return 1 389 390 output_lyr = CreateLayer(output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt) 391 if output_lyr is None: 392 return 1 393 else: 394 drv = output_ds.GetDriver() 395 396 if output_lyr_name is None: 397 cnt = output_ds.GetLayerCount() 398 if cnt != 1: 399 print('Result datasource has not a single layer, so you should specify its name with -output_lyr') 400 return 1 401 output_lyr = output_ds.GetLayer(0) 402 output_lyr_name = output_lyr.GetName() 403 else: 404 output_lyr = output_ds.GetLayerByName(output_lyr_name) 405 406 if output_lyr is None: 407 if EQUAL(drv.GetName(), "ESRI Shapefile") and \ 408 EQUAL(os.path.splitext(output_ds_name)[1], ".SHP") and \ 409 not EQUAL(output_lyr_name, os.path.splitext(os.path.basename(output_ds_name))[0]): 410 print('Cannot create layer "%s" in a shapefile called "%s"' % (output_lyr_name, output_ds_name)) 411 return 1 412 413 output_lyr = CreateLayer(output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt) 414 if output_lyr is None: 415 return 1 416 417 if overwrite: 418 cnt = output_ds.GetLayerCount() 419 iLayer = None # initialize in case there are no loop iterations 420 for iLayer in range(cnt): 421 poLayer = output_ds.GetLayer(iLayer) 422 if poLayer is not None \ 423 and poLayer.GetName() == output_lyr_name: 424 break 425 if iLayer != cnt: 426 if output_ds.DeleteLayer(iLayer) != 0: 427 print("DeleteLayer() failed when overwrite requested.") 428 return 1 429 430 output_lyr = CreateLayer(output_ds, output_lyr_name, srs, geom_type, lco, input_lyr, input_fields, method_lyr, method_fields, opt) 431 if output_lyr is None: 432 return 1 433 434 op = getattr(input_lyr, op_str) 435 if quiet_flag == 0: 436 ret = op(method_lyr, output_lyr, options=opt, callback=gdal.TermProgress_nocb) 437 else: 438 ret = op(method_lyr, output_lyr, options=opt) 439 440 input_ds = None 441 method_ds = None 442 output_ds = None 443 444 if ret != 0: 445 print('An error occurred during %s operation' % op_str) 446 return 1 447 448 return 0 449 450 451if __name__ == '__main__': 452 sys.exit(main(sys.argv)) 453