-
-
Save rouault/1c7b66b420d0ff665fbcb735e26e8664 to your computer and use it in GitHub Desktop.
osm2ogr.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# -*- coding: utf-8 -*- | |
#***************************************************************************** | |
# | |
# Project: OpenGIS Simple Features Reference Implementation | |
# Purpose: Python port of a simple client for translating between formats. | |
# Author: Even Rouault, <even dot rouault at mines dash paris dot org> | |
# | |
# Port from ogr2ogr.cpp whose author is Frank Warmerdam | |
# | |
#***************************************************************************** | |
# Copyright (c) 2010-2013, Even Rouault <even dot rouault at mines-paris dot org> | |
# Copyright (c) 1999, Frank Warmerdam | |
# | |
# 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 sys | |
import os | |
import stat | |
from osgeo import gdal | |
from osgeo import ogr | |
from osgeo import osr | |
############################################################################### | |
class ScaledProgressObject: | |
def __init__(self, min, max, cbk, cbk_data = None): | |
self.min = min | |
self.max = max | |
self.cbk = cbk | |
self.cbk_data = cbk_data | |
############################################################################### | |
def ScaledProgressFunc(pct, msg, data): | |
if data.cbk is None: | |
return True | |
return data.cbk(data.min + pct * (data.max - data.min), msg, data.cbk_data) | |
############################################################################### | |
def EQUAL(a, b): | |
return a.lower() == b.lower() | |
############################################################################### | |
# Redefinition of GDALTermProgress, so that autotest/pyscripts/test_osm2ogr_py.py | |
# can check that the progress bar is displayed | |
nLastTick = -1 | |
def TermProgress( dfComplete, pszMessage, pProgressArg ): | |
global nLastTick | |
nThisTick = (int) (dfComplete * 40.0) | |
if nThisTick < 0: | |
nThisTick = 0 | |
if nThisTick > 40: | |
nThisTick = 40 | |
# Have we started a new progress run? | |
if nThisTick < nLastTick and nLastTick >= 39: | |
nLastTick = -1 | |
if nThisTick <= nLastTick: | |
return True | |
while nThisTick > nLastTick: | |
nLastTick = nLastTick + 1 | |
if (nLastTick % 4) == 0: | |
sys.stdout.write('%d' % ((nLastTick / 4) * 10)) | |
else: | |
sys.stdout.write('.') | |
if nThisTick == 40: | |
print(" - done." ) | |
else: | |
sys.stdout.flush() | |
return True | |
class TargetLayerInfo: | |
def __init__(self): | |
self.poDstLayer = None | |
self.poCT = None | |
#self.papszTransformOptions = None | |
self.panMap = None | |
self.iSrcZField = None | |
class AssociatedLayers: | |
def __init__(self): | |
self.poSrcLayer = None | |
self.psInfo = None | |
#********************************************************************** | |
# main() | |
#********************************************************************** | |
bSkipFailures = False | |
nGroupTransactions = 200 | |
bPreserveFID = False | |
nFIDToFetch = ogr.NullFID | |
class Enum(set): | |
def __getattr__(self, name): | |
if name in self: | |
return name | |
raise AttributeError | |
GeomOperation = Enum(["NONE", "SEGMENTIZE", "SIMPLIFY_PRESERVE_TOPOLOGY"]) | |
def main(args = None, progress_func = TermProgress, progress_data = None): | |
global bSkipFailures | |
global nGroupTransactions | |
global bPreserveFID | |
global nFIDToFetch | |
pszFormat = "ESRI Shapefile" | |
pszDataSource = None | |
pszDestDataSource = None | |
papszLayers = [] | |
papszDSCO = [] | |
papszLCO = [] | |
bTransform = False | |
bAppend = False | |
bUpdate = False | |
bOverwrite = False | |
pszOutputSRSDef = None | |
pszSourceSRSDef = None | |
poOutputSRS = None | |
bNullifyOutputSRS = False | |
poSourceSRS = None | |
pszNewLayerName = None | |
pszWHERE = None | |
poSpatialFilter = None | |
pszSelect = None | |
papszSelFields = None | |
pszSQLStatement = None | |
eGType = -2 | |
bPromoteToMulti = False | |
eGeomOp = GeomOperation.NONE | |
dfGeomOpParam = 0 | |
papszFieldTypesToString = [] | |
bDisplayProgress = False | |
pfnProgress = None | |
pProgressArg = None | |
bClipSrc = False | |
bWrapDateline = False | |
poClipSrc = None | |
pszClipSrcDS = None | |
pszClipSrcSQL = None | |
pszClipSrcLayer = None | |
pszClipSrcWhere = None | |
poClipDst = None | |
pszClipDstDS = None | |
pszClipDstSQL = None | |
pszClipDstLayer = None | |
pszClipDstWhere = None | |
#pszSrcEncoding = None | |
#pszDstEncoding = None | |
bWrapDateline = False | |
bExplodeCollections = False | |
pszZField = None | |
nCoordDim = -1 | |
if args is None: | |
args = sys.argv | |
args = ogr.GeneralCmdLineProcessor( args ) | |
# -------------------------------------------------------------------- | |
# Processing command line arguments. | |
# -------------------------------------------------------------------- | |
if args is None: | |
return False | |
nArgc = len(args) | |
iArg = 1 | |
while iArg < nArgc: | |
if EQUAL(args[iArg],"-f") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszFormat = args[iArg] | |
elif EQUAL(args[iArg],"-dsco") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
papszDSCO.append(args[iArg] ) | |
elif EQUAL(args[iArg],"-lco") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
papszLCO.append(args[iArg] ) | |
elif EQUAL(args[iArg],"-preserve_fid"): | |
bPreserveFID = True | |
elif len(args[iArg]) >= 5 and EQUAL(args[iArg][0:5], "-skip"): | |
bSkipFailures = True | |
nGroupTransactions = 1 # #2409 | |
elif EQUAL(args[iArg],"-append"): | |
bAppend = True | |
bUpdate = True | |
elif EQUAL(args[iArg],"-overwrite"): | |
bOverwrite = True | |
bUpdate = True | |
elif EQUAL(args[iArg],"-update"): | |
bUpdate = True | |
elif EQUAL(args[iArg],"-fid") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
nFIDToFetch = int(args[iArg]) | |
elif EQUAL(args[iArg],"-sql") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszSQLStatement = args[iArg] | |
elif EQUAL(args[iArg],"-nln") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszNewLayerName = args[iArg] | |
elif EQUAL(args[iArg],"-nlt") and iArg < nArgc-1: | |
if EQUAL(args[iArg+1],"NONE"): | |
eGType = ogr.wkbNone | |
elif EQUAL(args[iArg+1],"GEOMETRY"): | |
eGType = ogr.wkbUnknown | |
elif EQUAL(args[iArg+1],"PROMOTE_TO_MULTI"): | |
bPromoteToMulti = True | |
elif EQUAL(args[iArg+1],"POINT"): | |
eGType = ogr.wkbPoint | |
elif EQUAL(args[iArg+1],"LINESTRING"): | |
eGType = ogr.wkbLineString | |
elif EQUAL(args[iArg+1],"POLYGON"): | |
eGType = ogr.wkbPolygon | |
elif EQUAL(args[iArg+1],"GEOMETRYCOLLECTION"): | |
eGType = ogr.wkbGeometryCollection | |
elif EQUAL(args[iArg+1],"MULTIPOINT"): | |
eGType = ogr.wkbMultiPoint | |
elif EQUAL(args[iArg+1],"MULTILINESTRING"): | |
eGType = ogr.wkbMultiLineString | |
elif EQUAL(args[iArg+1],"MULTIPOLYGON"): | |
eGType = ogr.wkbMultiPolygon | |
elif EQUAL(args[iArg+1],"GEOMETRY25D"): | |
eGType = ogr.wkbUnknown | ogr.wkb25DBit | |
elif EQUAL(args[iArg+1],"POINT25D"): | |
eGType = ogr.wkbPoint25D | |
elif EQUAL(args[iArg+1],"LINESTRING25D"): | |
eGType = ogr.wkbLineString25D | |
elif EQUAL(args[iArg+1],"POLYGON25D"): | |
eGType = ogr.wkbPolygon25D | |
elif EQUAL(args[iArg+1],"GEOMETRYCOLLECTION25D"): | |
eGType = ogr.wkbGeometryCollection25D | |
elif EQUAL(args[iArg+1],"MULTIPOINT25D"): | |
eGType = ogr.wkbMultiPoint25D | |
elif EQUAL(args[iArg+1],"MULTILINESTRING25D"): | |
eGType = ogr.wkbMultiLineString25D | |
elif EQUAL(args[iArg+1],"MULTIPOLYGON25D"): | |
eGType = ogr.wkbMultiPolygon25D | |
else: | |
print("-nlt %s: type not recognised." % args[iArg+1]) | |
return False | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-dim") and iArg < nArgc-1: | |
nCoordDim = int(args[iArg+1]) | |
if nCoordDim != 2 and nCoordDim != 3: | |
print("-dim %s: value not handled." % args[iArg+1]) | |
return False | |
iArg = iArg + 1 | |
elif (EQUAL(args[iArg],"-tg") or \ | |
EQUAL(args[iArg],"-gt")) and iArg < nArgc-1: | |
iArg = iArg + 1 | |
nGroupTransactions = int(args[iArg]) | |
elif EQUAL(args[iArg],"-s_srs") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszSourceSRSDef = args[iArg] | |
elif EQUAL(args[iArg],"-a_srs") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszOutputSRSDef = args[iArg] | |
if EQUAL(pszOutputSRSDef, "NULL") or \ | |
EQUAL(pszOutputSRSDef, "NONE"): | |
pszOutputSRSDef = None | |
bNullifyOutputSRS = True | |
elif EQUAL(args[iArg],"-t_srs") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszOutputSRSDef = args[iArg] | |
bTransform = True | |
elif EQUAL(args[iArg],"-spat") and iArg + 4 < nArgc: | |
oRing = ogr.Geometry(ogr.wkbLinearRing) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
poSpatialFilter = ogr.Geometry(ogr.wkbPolygon) | |
poSpatialFilter.AddGeometry(oRing) | |
iArg = iArg + 4 | |
elif EQUAL(args[iArg],"-where") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszWHERE = args[iArg] | |
elif EQUAL(args[iArg],"-select") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszSelect = args[iArg] | |
if pszSelect.find(',') != -1: | |
papszSelFields = pszSelect.split(',') | |
else: | |
papszSelFields = pszSelect.split(' ') | |
if papszSelFields[0] == '': | |
papszSelFields = [] | |
elif EQUAL(args[iArg],"-simplify") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
eGeomOp = GeomOperation.SIMPLIFY_PRESERVE_TOPOLOGY | |
dfGeomOpParam = float(args[iArg]) | |
elif EQUAL(args[iArg],"-segmentize") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
eGeomOp = GeomOperation.SEGMENTIZE | |
dfGeomOpParam = float(args[iArg]) | |
elif EQUAL(args[iArg],"-fieldTypeToString") and iArg < nArgc-1: | |
iArg = iArg + 1 | |
pszFieldTypeToString = args[iArg] | |
if pszFieldTypeToString.find(',') != -1: | |
tokens = pszFieldTypeToString.split(',') | |
else: | |
tokens = pszFieldTypeToString.split(' ') | |
for token in tokens: | |
if EQUAL(token,"Integer") or \ | |
EQUAL(token,"Real") or \ | |
EQUAL(token,"String") or \ | |
EQUAL(token,"Date") or \ | |
EQUAL(token,"Time") or \ | |
EQUAL(token,"DateTime") or \ | |
EQUAL(token,"Binary") or \ | |
EQUAL(token,"IntegerList") or \ | |
EQUAL(token,"RealList") or \ | |
EQUAL(token,"StringList"): | |
papszFieldTypesToString.append(token) | |
elif EQUAL(token,"All"): | |
papszFieldTypesToString = [ 'All' ] | |
break | |
else: | |
print("Unhandled type for fieldtypeasstring option : %s " % token) | |
return Usage() | |
elif EQUAL(args[iArg],"-progress"): | |
bDisplayProgress = True | |
#elif EQUAL(args[iArg],"-wrapdateline") ) | |
#{ | |
# bWrapDateline = True; | |
#} | |
# | |
elif EQUAL(args[iArg],"-clipsrc") and iArg < nArgc-1: | |
bClipSrc = True | |
if IsNumber(args[iArg+1]) and iArg < nArgc - 4: | |
oRing = ogr.Geometry(ogr.wkbLinearRing) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
poClipSrc = ogr.Geometry(ogr.wkbPolygon) | |
poClipSrc.AddGeometry(oRing) | |
iArg = iArg + 4 | |
elif (len(args[iArg+1]) >= 7 and EQUAL(args[iArg+1][0:7],"POLYGON") ) or \ | |
(len(args[iArg+1]) >= 12 and EQUAL(args[iArg+1][0:12],"MULTIPOLYGON") ) : | |
poClipSrc = ogr.CreateGeometryFromWkt(args[iArg+1]) | |
if poClipSrc is None: | |
print("FAILURE: Invalid geometry. Must be a valid POLYGON or MULTIPOLYGON WKT\n") | |
return Usage() | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg+1],"spat_extent"): | |
iArg = iArg + 1 | |
else: | |
pszClipSrcDS = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipsrcsql") and iArg < nArgc-1: | |
pszClipSrcSQL = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipsrclayer") and iArg < nArgc-1: | |
pszClipSrcLayer = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipsrcwhere") and iArg < nArgc-1: | |
pszClipSrcWhere = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipdst") and iArg < nArgc-1: | |
if IsNumber(args[iArg+1]) and iArg < nArgc - 4: | |
oRing = ogr.Geometry(ogr.wkbLinearRing) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+4]) ) | |
oRing.AddPoint_2D( float(args[iArg+3]), float(args[iArg+2]) ) | |
oRing.AddPoint_2D( float(args[iArg+1]), float(args[iArg+2]) ) | |
poClipDst = ogr.Geometry(ogr.wkbPolygon) | |
poClipDst.AddGeometry(oRing) | |
iArg = iArg + 4 | |
elif (len(args[iArg+1]) >= 7 and EQUAL(args[iArg+1][0:7],"POLYGON") ) or \ | |
(len(args[iArg+1]) >= 12 and EQUAL(args[iArg+1][0:12],"MULTIPOLYGON") ) : | |
poClipDst = ogr.CreateGeometryFromWkt(args[iArg+1]) | |
if poClipDst is None: | |
print("FAILURE: Invalid geometry. Must be a valid POLYGON or MULTIPOLYGON WKT\n") | |
return Usage() | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg+1],"spat_extent"): | |
iArg = iArg + 1 | |
else: | |
pszClipDstDS = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipdstsql") and iArg < nArgc-1: | |
pszClipDstSQL = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipdstlayer") and iArg < nArgc-1: | |
pszClipDstLayer = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-clipdstwhere") and iArg < nArgc-1: | |
pszClipDstWhere = args[iArg+1] | |
iArg = iArg + 1 | |
elif EQUAL(args[iArg],"-explodecollections"): | |
bExplodeCollections = True | |
elif EQUAL(args[iArg],"-zfield") and iArg < nArgc-1: | |
pszZField = args[iArg+1] | |
iArg = iArg + 1 | |
elif args[iArg][0] == '-': | |
return Usage() | |
elif pszDestDataSource is None: | |
pszDestDataSource = args[iArg] | |
elif pszDataSource is None: | |
pszDataSource = args[iArg] | |
else: | |
papszLayers.append (args[iArg] ) | |
iArg = iArg + 1 | |
if pszDataSource is None: | |
return Usage() | |
if bPreserveFID and bExplodeCollections: | |
print("FAILURE: cannot use -preserve_fid and -explodecollections at the same time\n\n") | |
return Usage() | |
if bClipSrc and pszClipSrcDS is not None: | |
poClipSrc = LoadGeometry(pszClipSrcDS, pszClipSrcSQL, pszClipSrcLayer, pszClipSrcWhere) | |
if poClipSrc is None: | |
print("FAILURE: cannot load source clip geometry\n" ) | |
return Usage() | |
elif bClipSrc and poClipSrc is None: | |
if poSpatialFilter is not None: | |
poClipSrc = poSpatialFilter.Clone() | |
if poClipSrc is None: | |
print("FAILURE: -clipsrc must be used with -spat option or a\n" + \ | |
"bounding box, WKT string or datasource must be specified\n") | |
return Usage() | |
if pszClipDstDS is not None: | |
poClipDst = LoadGeometry(pszClipDstDS, pszClipDstSQL, pszClipDstLayer, pszClipDstWhere) | |
if poClipDst is None: | |
print("FAILURE: cannot load dest clip geometry\n" ) | |
return Usage() | |
# -------------------------------------------------------------------- | |
# Open data source. | |
# -------------------------------------------------------------------- | |
poDS = ogr.Open( pszDataSource, False ) | |
# -------------------------------------------------------------------- | |
# Report failure | |
# -------------------------------------------------------------------- | |
if poDS is None: | |
print("FAILURE:\n" + \ | |
"Unable to open datasource `%s' with the following drivers." % pszDataSource) | |
for iDriver in range(ogr.GetDriverCount()): | |
print(" -> " + ogr.GetDriver(iDriver).GetName() ) | |
return False | |
# -------------------------------------------------------------------- | |
# Try opening the output datasource as an existing, writable | |
# -------------------------------------------------------------------- | |
poODS = None | |
poDriver = None | |
if bUpdate: | |
poODS = ogr.Open( pszDestDataSource, True ) | |
if poODS is None: | |
if bOverwrite or bAppend: | |
poODS = ogr.Open( pszDestDataSource, False ) | |
if poODS is None: | |
# the datasource doesn't exist at all | |
bUpdate = False | |
else: | |
poODS.delete() | |
poODS = None | |
if bUpdate: | |
print("FAILURE:\n" + | |
"Unable to open existing output datasource `%s'." % pszDestDataSource) | |
return False | |
elif len(papszDSCO) > 0: | |
print("WARNING: Datasource creation options ignored since an existing datasource\n" + \ | |
" being updated." ) | |
if poODS is not None: | |
poDriver = poODS.GetDriver() | |
# -------------------------------------------------------------------- | |
# Find the output driver. | |
# -------------------------------------------------------------------- | |
if not bUpdate: | |
poDriver = ogr.GetDriverByName(pszFormat) | |
if poDriver is None: | |
print("Unable to find driver `%s'." % pszFormat) | |
print( "The following drivers are available:" ) | |
for iDriver in range(ogr.GetDriverCount()): | |
print(" -> %s" % ogr.GetDriver(iDriver).GetName() ) | |
return False | |
if poDriver.TestCapability( ogr.ODrCCreateDataSource ) == False: | |
print( "%s driver does not support data source creation." % pszFormat) | |
return False | |
# -------------------------------------------------------------------- | |
# Special case to improve user experience when translating | |
# a datasource with multiple layers into a shapefile. If the | |
# user gives a target datasource with .shp and it does not exist, | |
# the shapefile driver will try to create a file, but this is not | |
# appropriate because here we have several layers, so create | |
# a directory instead. | |
# -------------------------------------------------------------------- | |
if EQUAL(poDriver.GetName(), "ESRI Shapefile") and \ | |
pszSQLStatement is None and \ | |
(len(papszLayers) > 1 or \ | |
(len(papszLayers) == 0 and poDS.GetLayerCount() > 1)) and \ | |
pszNewLayerName is None and \ | |
EQUAL(os.path.splitext(pszDestDataSource)[1], ".SHP") : | |
try: | |
os.stat(pszDestDataSource) | |
except: | |
try: | |
# decimal 493 = octal 0755. Python 3 needs 0o755, but | |
# this syntax is only supported by Python >= 2.6 | |
os.mkdir(pszDestDataSource, 493) | |
except: | |
print("Failed to create directory %s\n" | |
"for shapefile datastore.\n" % pszDestDataSource ) | |
return False | |
# -------------------------------------------------------------------- | |
# Create the output data source. | |
# -------------------------------------------------------------------- | |
poODS = poDriver.CreateDataSource( pszDestDataSource, options = papszDSCO ) | |
if poODS is None: | |
print( "%s driver failed to create %s" % (pszFormat, pszDestDataSource )) | |
return False | |
# -------------------------------------------------------------------- | |
# Parse the output SRS definition if possible. | |
# -------------------------------------------------------------------- | |
if pszOutputSRSDef is not None: | |
poOutputSRS = osr.SpatialReference() | |
if poOutputSRS.SetFromUserInput( pszOutputSRSDef ) != 0: | |
print( "Failed to process SRS definition: %s" % pszOutputSRSDef ) | |
return False | |
# -------------------------------------------------------------------- | |
# Parse the source SRS definition if possible. | |
# -------------------------------------------------------------------- | |
if pszSourceSRSDef is not None: | |
poSourceSRS = osr.SpatialReference() | |
if poSourceSRS.SetFromUserInput( pszSourceSRSDef ) != 0: | |
print( "Failed to process SRS definition: %s" % pszSourceSRSDef ) | |
return False | |
# -------------------------------------------------------------------- | |
# For OSM file. | |
# -------------------------------------------------------------------- | |
bSrcIsOSM = poDS.GetDriver() is not None and \ | |
poDS.GetDriver().GetName() == "OSM" | |
nSrcFileSize = 0 | |
if bSrcIsOSM and poDS.GetName() != "/vsistdin/": | |
sStat = gdal.VSIStatL(poDS.GetName()) | |
if sStat is not None: | |
nSrcFileSize = sStat.size | |
# -------------------------------------------------------------------- | |
# Special case for -sql clause. No source layers required. | |
# -------------------------------------------------------------------- | |
if pszSQLStatement is not None: | |
if pszWHERE is not None: | |
print( "-where clause ignored in combination with -sql." ) | |
if len(papszLayers) > 0: | |
print( "layer names ignored in combination with -sql." ) | |
poResultSet = poDS.ExecuteSQL( pszSQLStatement, poSpatialFilter, \ | |
None ) | |
if poResultSet is not None: | |
nCountLayerFeatures = 0 | |
if bDisplayProgress: | |
if bSrcIsOSM: | |
pfnProgress = progress_func | |
pProgressArg = progress_data | |
elif not poResultSet.TestCapability(ogr.OLCFastFeatureCount): | |
print( "Progress turned off as fast feature count is not available.") | |
bDisplayProgress = False | |
else: | |
nCountLayerFeatures = poResultSet.GetFeatureCount() | |
pfnProgress = progress_func | |
pProgressArg = progress_data | |
# -------------------------------------------------------------------- | |
# Special case to improve user experience when translating into | |
# single file shapefile and source has only one layer, and that | |
# the layer name isn't specified | |
# -------------------------------------------------------------------- | |
if EQUAL(poDriver.GetName(), "ESRI Shapefile") and \ | |
pszNewLayerName is None: | |
try: | |
mode = os.stat(pszDestDataSource).st_mode | |
if (mode & stat.S_IFDIR) == 0: | |
pszNewLayerName = os.path.splitext(os.path.basename(pszDestDataSource))[0] | |
except: | |
pass | |
psInfo = SetupTargetLayer( poDS, \ | |
poResultSet, | |
poODS, \ | |
papszLCO, \ | |
pszNewLayerName, \ | |
bTransform, \ | |
poOutputSRS, \ | |
bNullifyOutputSRS, \ | |
poSourceSRS, \ | |
papszSelFields, \ | |
bAppend, eGType, bPromoteToMulti, nCoordDim, bOverwrite, \ | |
papszFieldTypesToString, \ | |
bWrapDateline, \ | |
bExplodeCollections, \ | |
pszZField, \ | |
pszWHERE ) | |
poResultSet.ResetReading() | |
if psInfo is None or not TranslateLayer( psInfo, poDS, poResultSet, poODS, \ | |
poOutputSRS, bNullifyOutputSRS, \ | |
eGType, bPromoteToMulti, nCoordDim, \ | |
eGeomOp, dfGeomOpParam, \ | |
nCountLayerFeatures, \ | |
poClipSrc, poClipDst, \ | |
bExplodeCollections, \ | |
nSrcFileSize, None, \ | |
pfnProgress, pProgressArg ): | |
print( | |
"Terminating translation prematurely after failed\n" + \ | |
"translation from sql statement." ) | |
return False | |
poDS.ReleaseResultSet( poResultSet ) | |
# -------------------------------------------------------------------- | |
# Special case for layer interleaving mode. | |
# -------------------------------------------------------------------- | |
elif bSrcIsOSM and gdal.GetConfigOption("OGR_INTERLEAVED_READING", None) is None: | |
gdal.SetConfigOption("OGR_INTERLEAVED_READING", "YES") | |
#if (bSplitListFields) | |
#{ | |
# fprintf( stderr, "FAILURE: -splitlistfields not supported in this mode\n" ); | |
# exit( 1 ); | |
#} | |
nSrcLayerCount = poDS.GetLayerCount() | |
pasAssocLayers = [ AssociatedLayers() for i in range(nSrcLayerCount) ] | |
# -------------------------------------------------------------------- | |
# Special case to improve user experience when translating into | |
# single file shapefile and source has only one layer, and that | |
# the layer name isn't specified | |
# -------------------------------------------------------------------- | |
if EQUAL(poDriver.GetName(), "ESRI Shapefile") and \ | |
(len(papszLayers) == 1 or nSrcLayerCount == 1) and pszNewLayerName is None: | |
try: | |
mode = os.stat(pszDestDataSource).st_mode | |
if (mode & stat.S_IFDIR) == 0: | |
pszNewLayerName = os.path.splitext(os.path.basename(pszDestDataSource))[0] | |
except: | |
pass | |
if bDisplayProgress and bSrcIsOSM: | |
pfnProgress = progress_func | |
pProgressArg = progress_data | |
# -------------------------------------------------------------------- | |
# If no target layer specified, use all source layers. | |
# -------------------------------------------------------------------- | |
if len(papszLayers) == 0: | |
papszLayers = [ None for i in range(nSrcLayerCount) ] | |
for iLayer in range(nSrcLayerCount): | |
poLayer = poDS.GetLayer(iLayer) | |
if poLayer is None: | |
print("FAILURE: Couldn't fetch advertised layer %d!" % iLayer) | |
return False | |
papszLayers[iLayer] = poLayer.GetName() | |
else: | |
if bSrcIsOSM: | |
osInterestLayers = "SET interest_layers =" | |
for iLayer in range(len(papszLayers)): | |
if iLayer != 0: | |
osInterestLayers = osInterestLayers + "," | |
osInterestLayers = osInterestLayers + papszLayers[iLayer] | |
poDS.ExecuteSQL(osInterestLayers, None, None) | |
# -------------------------------------------------------------------- | |
# First pass to set filters and create target layers. | |
# -------------------------------------------------------------------- | |
for iLayer in range(nSrcLayerCount): | |
poLayer = poDS.GetLayer(iLayer) | |
if poLayer is None: | |
print("FAILURE: Couldn't fetch advertised layer %d!" % iLayer) | |
return False | |
pasAssocLayers[iLayer].poSrcLayer = poLayer | |
if CSLFindString(papszLayers, poLayer.GetName()) >= 0: | |
if pszWHERE is not None: | |
if poLayer.SetAttributeFilter( pszWHERE ) != 0: | |
print("FAILURE: SetAttributeFilter(%s) on layer '%s' failed.\n" % (pszWHERE, poLayer.GetName()) ) | |
if not bSkipFailures: | |
return False | |
if poSpatialFilter is not None: | |
poLayer.SetSpatialFilter( poSpatialFilter ) | |
psInfo = SetupTargetLayer( poDS, \ | |
poLayer, \ | |
poODS, \ | |
papszLCO, \ | |
pszNewLayerName, \ | |
bTransform, \ | |
poOutputSRS, \ | |
bNullifyOutputSRS, \ | |
poSourceSRS, \ | |
papszSelFields, \ | |
bAppend, eGType, bPromoteToMulti, nCoordDim, bOverwrite, \ | |
papszFieldTypesToString, \ | |
bWrapDateline, \ | |
bExplodeCollections, \ | |
pszZField, \ | |
pszWHERE ) | |
if psInfo is None and not bSkipFailures: | |
return False | |
pasAssocLayers[iLayer].psInfo = psInfo | |
else: | |
pasAssocLayers[iLayer].psInfo = None | |
# -------------------------------------------------------------------- | |
# Second pass to process features in a interleaved layer mode. | |
# -------------------------------------------------------------------- | |
bHasLayersNonEmpty = True | |
while bHasLayersNonEmpty: | |
bHasLayersNonEmpty = False | |
for iLayer in range(nSrcLayerCount): | |
poLayer = pasAssocLayers[iLayer].poSrcLayer | |
psInfo = pasAssocLayers[iLayer].psInfo | |
anReadFeatureCount = [0] | |
if psInfo is not None: | |
if not TranslateLayer(psInfo, poDS, poLayer, poODS, \ | |
poOutputSRS, bNullifyOutputSRS, \ | |
eGType, bPromoteToMulti, nCoordDim, \ | |
eGeomOp, dfGeomOpParam, \ | |
0, \ | |
poClipSrc, poClipDst, \ | |
bExplodeCollections, \ | |
nSrcFileSize, \ | |
anReadFeatureCount, \ | |
pfnProgress, pProgressArg ) \ | |
and not bSkipFailures: | |
print( | |
"Terminating translation prematurely after failed\n" + \ | |
"translation of layer " + poLayer.GetName() + " (use -skipfailures to skip errors)") | |
return False | |
else: | |
# No matching target layer : just consumes the features | |
poFeature = poLayer.GetNextFeature() | |
while poFeature is not None: | |
anReadFeatureCount[0] = anReadFeatureCount[0] + 1 | |
poFeature = poLayer.GetNextFeature() | |
if anReadFeatureCount[0] != 0: | |
bHasLayersNonEmpty = True | |
else: | |
nLayerCount = 0 | |
papoLayers = [] | |
# -------------------------------------------------------------------- | |
# Process each data source layer. | |
# -------------------------------------------------------------------- | |
if len(papszLayers) == 0: | |
nLayerCount = poDS.GetLayerCount() | |
papoLayers = [None for i in range(nLayerCount)] | |
iLayer = 0 | |
for iLayer in range(nLayerCount): | |
poLayer = poDS.GetLayer(iLayer) | |
if poLayer is None: | |
print("FAILURE: Couldn't fetch advertised layer %d!" % iLayer) | |
return False | |
papoLayers[iLayer] = poLayer | |
iLayer = iLayer + 1 | |
# -------------------------------------------------------------------- | |
# Process specified data source layers. | |
# -------------------------------------------------------------------- | |
else: | |
nLayerCount = len(papszLayers) | |
papoLayers = [None for i in range(nLayerCount)] | |
iLayer = 0 | |
for layername in papszLayers: | |
poLayer = poDS.GetLayerByName(layername) | |
if poLayer is None: | |
print("FAILURE: Couldn't fetch advertised layer %s!" % layername) | |
return False | |
papoLayers[iLayer] = poLayer | |
iLayer = iLayer + 1 | |
panLayerCountFeatures = [0 for i in range(nLayerCount)] | |
nCountLayersFeatures = 0 | |
nAccCountFeatures = 0 | |
# First pass to apply filters and count all features if necessary | |
for iLayer in range(nLayerCount): | |
poLayer = papoLayers[iLayer] | |
if pszWHERE is not None: | |
if poLayer.SetAttributeFilter( pszWHERE ) != 0: | |
print("FAILURE: SetAttributeFilter(%s) failed." % pszWHERE) | |
if not bSkipFailures: | |
return False | |
if poSpatialFilter is not None: | |
poLayer.SetSpatialFilter( poSpatialFilter ) | |
if bDisplayProgress and not bSrcIsOSM: | |
if not poLayer.TestCapability(ogr.OLCFastFeatureCount): | |
print("Progress turned off as fast feature count is not available.") | |
bDisplayProgress = False | |
else: | |
panLayerCountFeatures[iLayer] = poLayer.GetFeatureCount() | |
nCountLayersFeatures += panLayerCountFeatures[iLayer] | |
# Second pass to do the real job | |
for iLayer in range(nLayerCount): | |
poLayer = papoLayers[iLayer] | |
if bDisplayProgress: | |
if bSrcIsOSM: | |
pfnProgress = progress_func | |
pProgressArg = progress_data | |
else: | |
pfnProgress = ScaledProgressFunc | |
pProgressArg = ScaledProgressObject( \ | |
nAccCountFeatures * 1.0 / nCountLayersFeatures, \ | |
(nAccCountFeatures + panLayerCountFeatures[iLayer]) * 1.0 / nCountLayersFeatures, \ | |
progress_func, progress_data) | |
nAccCountFeatures += panLayerCountFeatures[iLayer] | |
# -------------------------------------------------------------------- | |
# Special case to improve user experience when translating into | |
# single file shapefile and source has only one layer, and that | |
# the layer name isn't specified | |
# -------------------------------------------------------------------- | |
if EQUAL(poDriver.GetName(), "ESRI Shapefile") and \ | |
nLayerCount == 1 and pszNewLayerName is None: | |
try: | |
mode = os.stat(pszDestDataSource).st_mode | |
if (mode & stat.S_IFDIR) == 0: | |
pszNewLayerName = os.path.splitext(os.path.basename(pszDestDataSource))[0] | |
except: | |
pass | |
psInfo = SetupTargetLayer( poDS, \ | |
poLayer, \ | |
poODS, \ | |
papszLCO, \ | |
pszNewLayerName, \ | |
bTransform, \ | |
poOutputSRS, \ | |
bNullifyOutputSRS, \ | |
poSourceSRS, \ | |
papszSelFields, \ | |
bAppend, eGType, bPromoteToMulti, nCoordDim, bOverwrite, \ | |
papszFieldTypesToString, \ | |
bWrapDateline, \ | |
bExplodeCollections, \ | |
pszZField, \ | |
pszWHERE ) | |
poLayer.ResetReading() | |
if (psInfo is None or \ | |
not TranslateLayer( psInfo, poDS, poLayer, poODS, \ | |
poOutputSRS, bNullifyOutputSRS, \ | |
eGType, bPromoteToMulti, nCoordDim, \ | |
eGeomOp, dfGeomOpParam, \ | |
panLayerCountFeatures[iLayer], \ | |
poClipSrc, poClipDst, \ | |
bExplodeCollections, \ | |
nSrcFileSize, None, \ | |
pfnProgress, pProgressArg )) \ | |
and not bSkipFailures: | |
print( | |
"Terminating translation prematurely after failed\n" + \ | |
"translation of layer " + poLayer.GetLayerDefn().GetName() + " (use -skipfailures to skip errors)") | |
return False | |
# -------------------------------------------------------------------- | |
# Close down. | |
# -------------------------------------------------------------------- | |
# We must explicitly destroy the output dataset in order the file | |
# to be properly closed ! | |
poODS.Destroy() | |
poDS.Destroy() | |
return True | |
#********************************************************************** | |
# Usage() | |
#********************************************************************** | |
def Usage(): | |
print( "Usage: osm2ogr [--help-general] [-skipfailures] [-append] [-update] [-gt n]\n" + \ | |
" [-select field_list] [-where restricted_where] \n" + \ | |
" [-progress] [-sql <sql statement>] \n" + \ | |
" [-spat xmin ymin xmax ymax] [-preserve_fid] [-fid FID]\n" + \ | |
" [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]\n" + \ | |
" [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]\n" + \ | |
" [-simplify tolerance]\n" + \ | |
#// " [-segmentize max_dist] [-fieldTypeToString All|(type1[,type2]*)]\n" + \ | |
" [-fieldTypeToString All|(type1[,type2]*)] [-explodecollections] \n" + \ | |
" dst_datasource_name src_datasource_name\n" + \ | |
" [-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3] [layer [layer ...]]\n" + \ | |
"\n" + \ | |
" -f format_name: output file format name, possible values are:") | |
for iDriver in range(ogr.GetDriverCount()): | |
poDriver = ogr.GetDriver(iDriver) | |
if poDriver.TestCapability( ogr.ODrCCreateDataSource ): | |
print( " -f \"" + poDriver.GetName() + "\"" ) | |
print( " -append: Append to existing layer instead of creating new if it exists\n" + \ | |
" -overwrite: delete the output layer and recreate it empty\n" + \ | |
" -update: Open existing output datasource in update mode\n" + \ | |
" -progress: Display progress on terminal. Only works if input layers have the \"fast feature count\" capability\n" + \ | |
" -select field_list: Comma-delimited list of fields from input layer to\n" + \ | |
" copy to the new layer (defaults to all)\n" + \ | |
" -where restricted_where: Attribute query (like SQL WHERE)\n" + \ | |
" -sql statement: Execute given SQL statement and save result.\n" + \ | |
" -skipfailures: skip features or layers that fail to convert\n" + \ | |
" -gt n: group n features per transaction (default 200)\n" + \ | |
" -spat xmin ymin xmax ymax: spatial query extents\n" + \ | |
" -simplify tolerance: distance tolerance for simplification.\n" + \ | |
#//" -segmentize max_dist: maximum distance between 2 nodes.\n" + \ | |
#//" Used to create intermediate points\n" + \ | |
" -dsco NAME=VALUE: Dataset creation option (format specific)\n" + \ | |
" -lco NAME=VALUE: Layer creation option (format specific)\n" + \ | |
" -nln name: Assign an alternate name to the new layer\n" + \ | |
" -nlt type: Force a geometry type for new layer. One of NONE, GEOMETRY,\n" + \ | |
" POINT, LINESTRING, POLYGON, GEOMETRYCOLLECTION, MULTIPOINT,\n" + \ | |
" MULTIPOLYGON, or MULTILINESTRING. Add \"25D\" for 3D layers.\n" + \ | |
" Default is type of source layer.\n" + \ | |
" -dim dimension: Force the coordinate dimension to the specified value.\n" + \ | |
" -fieldTypeToString type1,...: Converts fields of specified types to\n" + \ | |
" fields of type string in the new layer. Valid types are : \n" + \ | |
" Integer, Real, String, Date, Time, DateTime, Binary, IntegerList, RealList,\n" + \ | |
" StringList. Special value All can be used to convert all fields to strings.") | |
print(" -a_srs srs_def: Assign an output SRS\n" | |
" -t_srs srs_def: Reproject/transform to this SRS on output\n" | |
" -s_srs srs_def: Override source SRS\n" | |
"\n" | |
" Srs_def can be a full WKT definition (hard to escape properly),\n" | |
" or a well known definition (i.e. EPSG:4326) or a file with a WKT\n" | |
" definition." ) | |
return False | |
def CSLFindString(v, mystr): | |
i = 0 | |
for strIter in v: | |
if EQUAL(strIter, mystr): | |
return i | |
i = i + 1 | |
return -1 | |
def IsNumber( pszStr): | |
try: | |
(float)(pszStr) | |
return True | |
except: | |
return False | |
def LoadGeometry( pszDS, pszSQL, pszLyr, pszWhere): | |
poGeom = None | |
poDS = ogr.Open( pszDS, False ) | |
if poDS is None: | |
return None | |
if pszSQL is not None: | |
poLyr = poDS.ExecuteSQL( pszSQL, None, None ) | |
elif pszLyr is not None: | |
poLyr = poDS.GetLayerByName(pszLyr) | |
else: | |
poLyr = poDS.GetLayer(0) | |
if poLyr is None: | |
print("Failed to identify source layer from datasource.") | |
poDS.Destroy() | |
return None | |
if pszWhere is not None: | |
poLyr.SetAttributeFilter(pszWhere) | |
poFeat = poLyr.GetNextFeature() | |
while poFeat is not None: | |
poSrcGeom = poFeat.GetGeometryRef() | |
if poSrcGeom is not None: | |
eType = wkbFlatten(poSrcGeom.GetGeometryType()) | |
if poGeom is None: | |
poGeom = ogr.Geometry( ogr.wkbMultiPolygon ) | |
if eType == ogr.wkbPolygon: | |
poGeom.AddGeometry( poSrcGeom ) | |
elif eType == ogr.wkbMultiPolygon: | |
for iGeom in range(poSrcGeom.GetGeometryCount()): | |
poGeom.AddGeometry(poSrcGeom.GetGeometryRef(iGeom) ) | |
else: | |
print("ERROR: Geometry not of polygon type." ) | |
if pszSQL is not None: | |
poDS.ReleaseResultSet( poLyr ) | |
poDS.Destroy() | |
return None | |
poFeat = poLyr.GetNextFeature() | |
if pszSQL is not None: | |
poDS.ReleaseResultSet( poLyr ) | |
poDS.Destroy() | |
return poGeom | |
def wkbFlatten(x): | |
return x & (~ogr.wkb25DBit) | |
#********************************************************************** | |
# SetZ() | |
#********************************************************************** | |
def SetZ (poGeom, dfZ ): | |
if poGeom is None: | |
return | |
eGType = wkbFlatten(poGeom.GetGeometryType()) | |
if eGType == ogr.wkbPoint: | |
poGeom.SetPoint(0, poGeom.GetX(), poGeom.GetY(), dfZ) | |
elif eGType == ogr.wkbLineString or \ | |
eGType == ogr.wkbLinearRing: | |
for i in range(poGeom.GetPointCount()): | |
poGeom.SetPoint(i, poGeom.GetX(i), poGeom.GetY(i), dfZ) | |
elif eGType == ogr.wkbPolygon or \ | |
eGType == ogr.wkbMultiPoint or \ | |
eGType == ogr.wkbMultiLineString or \ | |
eGType == ogr.wkbMultiPolygon or \ | |
eGType == ogr.wkbGeometryCollection: | |
for i in range(poGeom.GetGeometryCount()): | |
SetZ(poGeom.GetGeometryRef(i), dfZ) | |
#********************************************************************** | |
# SetupTargetLayer() | |
#********************************************************************** | |
def SetupTargetLayer( poSrcDS, poSrcLayer, poDstDS, papszLCO, pszNewLayerName, \ | |
bTransform, poOutputSRS, bNullifyOutputSRS, poSourceSRS, papszSelFields, \ | |
bAppend, eGType, bPromoteToMulti, nCoordDim, bOverwrite, \ | |
papszFieldTypesToString, bWrapDateline, \ | |
bExplodeCollections, pszZField, pszWHERE) : | |
if pszNewLayerName is None: | |
pszNewLayerName = poSrcLayer.GetLayerDefn().GetName() | |
# -------------------------------------------------------------------- | |
# Setup coordinate transformation if we need it. | |
# -------------------------------------------------------------------- | |
poCT = None | |
if bTransform: | |
if poSourceSRS is None: | |
poSourceSRS = poSrcLayer.GetSpatialRef() | |
if poSourceSRS is None: | |
print("Can't transform coordinates, source layer has no\n" + \ | |
"coordinate system. Use -s_srs to set one." ) | |
return None | |
poCT = osr.CoordinateTransformation( poSourceSRS, poOutputSRS ) | |
if gdal.GetLastErrorMsg().find( 'Unable to load PROJ.4 library' ) != -1: | |
poCT = None | |
if poCT is None: | |
pszWKT = None | |
print("Failed to create coordinate transformation between the\n" + \ | |
"following coordinate systems. This may be because they\n" + \ | |
"are not transformable, or because projection services\n" + \ | |
"(PROJ.4 DLL/.so) could not be loaded." ) | |
pszWKT = poSourceSRS.ExportToPrettyWkt( 0 ) | |
print( "Source:\n" + pszWKT ) | |
pszWKT = poOutputSRS.ExportToPrettyWkt( 0 ) | |
print( "Target:\n" + pszWKT ) | |
return None | |
# -------------------------------------------------------------------- | |
# Get other info. | |
# -------------------------------------------------------------------- | |
poSrcFDefn = poSrcLayer.GetLayerDefn() | |
if poOutputSRS is None and not bNullifyOutputSRS: | |
poOutputSRS = poSrcLayer.GetSpatialRef() | |
# -------------------------------------------------------------------- | |
# Find the layer. | |
# -------------------------------------------------------------------- | |
# GetLayerByName() can instantiate layers that would have been | |
# 'hidden' otherwise, for example, non-spatial tables in a | |
# PostGIS-enabled database, so this apparently useless command is | |
# not useless. (#4012) | |
gdal.PushErrorHandler('CPLQuietErrorHandler') | |
poDstLayer = poDstDS.GetLayerByName(pszNewLayerName) | |
gdal.PopErrorHandler() | |
gdal.ErrorReset() | |
iLayer = -1 | |
if poDstLayer is not None: | |
nLayerCount = poDstDS.GetLayerCount() | |
for iLayer in range(nLayerCount): | |
poLayer = poDstDS.GetLayer(iLayer) | |
# The .cpp version compares on pointers directly, but we cannot | |
# do this with swig object, so just compare the names. | |
if poLayer is not None \ | |
and poLayer.GetName() == poDstLayer.GetName(): | |
break | |
if (iLayer == nLayerCount): | |
# Shouldn't happen with an ideal driver | |
poDstLayer = None | |
# -------------------------------------------------------------------- | |
# If the user requested overwrite, and we have the layer in | |
# question we need to delete it now so it will get recreated | |
# (overwritten). | |
# -------------------------------------------------------------------- | |
if poDstLayer is not None and bOverwrite: | |
if poDstDS.DeleteLayer( iLayer ) != 0: | |
print("DeleteLayer() failed when overwrite requested." ) | |
return None | |
poDstLayer = None | |
# -------------------------------------------------------------------- | |
# If the layer does not exist, then create it. | |
# -------------------------------------------------------------------- | |
if poDstLayer is None: | |
if eGType == -2: | |
eGType = poSrcFDefn.GetGeomType() | |
n25DBit = eGType & ogr.wkb25DBit | |
if bPromoteToMulti: | |
if wkbFlatten(eGType) == ogr.wkbLineString: | |
eGType = ogr.wkbMultiLineString | n25DBit | |
elif wkbFlatten(eGType) == ogr.wkbPolygon: | |
eGType = ogr.wkbMultiPolygon | n25DBit | |
if bExplodeCollections: | |
if wkbFlatten(eGType) == ogr.wkbMultiPoint: | |
eGType = ogr.wkbPoint | n25DBit | |
elif wkbFlatten(eGType) == ogr.wkbMultiLineString: | |
eGType = ogr.wkbLineString | n25DBit | |
elif wkbFlatten(eGType) == ogr.wkbMultiPolygon: | |
eGType = ogr.wkbPolygon | n25DBit | |
elif wkbFlatten(eGType) == ogr.wkbGeometryCollection: | |
eGType = ogr.wkbUnknown | n25DBit | |
if pszZField is not None: | |
eGType = eGType | ogr.wkb25DBit | |
if nCoordDim == 2: | |
eGType = eGType & ~ogr.wkb25DBit | |
elif nCoordDim == 3: | |
eGType = eGType | ogr.wkb25DBit | |
if poDstDS.TestCapability( ogr.ODsCCreateLayer ) == False: | |
print("Layer " + pszNewLayerName + "not found, and CreateLayer not supported by driver.") | |
return None | |
gdal.ErrorReset() | |
poDstLayer = poDstDS.CreateLayer( pszNewLayerName, poOutputSRS, \ | |
eGType, papszLCO ) | |
if poDstLayer is None: | |
return None | |
bAppend = False | |
# -------------------------------------------------------------------- | |
# Otherwise we will append to it, if append was requested. | |
# -------------------------------------------------------------------- | |
elif not bAppend: | |
print("FAILED: Layer " + pszNewLayerName + "already exists, and -append not specified.\n" + \ | |
" Consider using -append, or -overwrite.") | |
return None | |
else: | |
if len(papszLCO) > 0: | |
print("WARNING: Layer creation options ignored since an existing layer is\n" + \ | |
" being appended to." ) | |
# -------------------------------------------------------------------- | |
# Add fields. Default to copy all field. | |
# If only a subset of all fields requested, then output only | |
# the selected fields, and in the order that they were | |
# selected. | |
# -------------------------------------------------------------------- | |
# Initialize the index-to-index map to -1's | |
nSrcFieldCount = poSrcFDefn.GetFieldCount() | |
panMap = [ -1 for i in range(nSrcFieldCount) ] | |
poDstFDefn = poDstLayer.GetLayerDefn() | |
if papszSelFields is not None and not bAppend: | |
nDstFieldCount = 0 | |
if poDstFDefn is not None: | |
nDstFieldCount = poDstFDefn.GetFieldCount() | |
for iField in range(len(papszSelFields)): | |
iSrcField = poSrcFDefn.GetFieldIndex(papszSelFields[iField]) | |
if iSrcField >= 0: | |
poSrcFieldDefn = poSrcFDefn.GetFieldDefn(iSrcField) | |
oFieldDefn = ogr.FieldDefn( poSrcFieldDefn.GetNameRef(), | |
poSrcFieldDefn.GetType() ) | |
oFieldDefn.SetWidth( poSrcFieldDefn.GetWidth() ) | |
oFieldDefn.SetPrecision( poSrcFieldDefn.GetPrecision() ) | |
if papszFieldTypesToString is not None and \ | |
(CSLFindString(papszFieldTypesToString, "All") != -1 or \ | |
CSLFindString(papszFieldTypesToString, \ | |
ogr.GetFieldTypeName(poSrcFieldDefn.GetType())) != -1): | |
oFieldDefn.SetType(ogr.OFTString) | |
# The field may have been already created at layer creation | |
iDstField = -1 | |
if poDstFDefn is not None: | |
iDstField = poDstFDefn.GetFieldIndex(oFieldDefn.GetNameRef()) | |
if iDstField >= 0: | |
panMap[iSrcField] = iDstField | |
elif poDstLayer.CreateField( oFieldDefn ) == 0: | |
# now that we've created a field, GetLayerDefn() won't return NULL | |
if poDstFDefn is None: | |
poDstFDefn = poDstLayer.GetLayerDefn() | |
# Sanity check : if it fails, the driver is buggy | |
if poDstFDefn is not None and \ | |
poDstFDefn.GetFieldCount() != nDstFieldCount + 1: | |
print("The output driver has claimed to have added the %s field, but it did not!" % oFieldDefn.GetNameRef() ) | |
else: | |
panMap[iSrcField] = nDstFieldCount | |
nDstFieldCount = nDstFieldCount + 1 | |
else: | |
print("Field '" + papszSelFields[iField] + "' not found in source layer.") | |
if not bSkipFailures: | |
return None | |
# -------------------------------------------------------------------- | |
# Use SetIgnoredFields() on source layer if available | |
# -------------------------------------------------------------------- | |
# Here we differ from the osm2ogr.cpp implementation since the OGRFeatureQuery | |
# isn't mapped to swig. So in that case just don't use SetIgnoredFields() | |
# to avoid issue raised in #4015 | |
if poSrcLayer.TestCapability(ogr.OLCIgnoreFields) and pszWHERE is None: | |
papszIgnoredFields = [] | |
for iSrcField in range(nSrcFieldCount): | |
pszFieldName = poSrcFDefn.GetFieldDefn(iSrcField).GetNameRef() | |
bFieldRequested = False | |
for iField in range(len(papszSelFields)): | |
if EQUAL(pszFieldName, papszSelFields[iField]): | |
bFieldRequested = True | |
break | |
if pszZField is not None and EQUAL(pszFieldName, pszZField): | |
bFieldRequested = True | |
# If source field not requested, add it to ignored files list | |
if not bFieldRequested: | |
papszIgnoredFields.append(pszFieldName) | |
poSrcLayer.SetIgnoredFields(papszIgnoredFields) | |
elif not bAppend: | |
nDstFieldCount = 0 | |
if poDstFDefn is not None: | |
nDstFieldCount = poDstFDefn.GetFieldCount() | |
for iField in range(nSrcFieldCount): | |
poSrcFieldDefn = poSrcFDefn.GetFieldDefn(iField) | |
if poSrcFieldDefn.GetNameRef() == 'other_tags': | |
continue | |
oFieldDefn = ogr.FieldDefn( poSrcFieldDefn.GetNameRef(), | |
poSrcFieldDefn.GetType() ) | |
oFieldDefn.SetWidth( poSrcFieldDefn.GetWidth() ) | |
oFieldDefn.SetPrecision( poSrcFieldDefn.GetPrecision() ) | |
if papszFieldTypesToString is not None and \ | |
(CSLFindString(papszFieldTypesToString, "All") != -1 or \ | |
CSLFindString(papszFieldTypesToString, \ | |
ogr.GetFieldTypeName(poSrcFieldDefn.GetType())) != -1): | |
oFieldDefn.SetType(ogr.OFTString) | |
# The field may have been already created at layer creation | |
iDstField = -1 | |
if poDstFDefn is not None: | |
iDstField = poDstFDefn.GetFieldIndex(oFieldDefn.GetNameRef()) | |
if iDstField >= 0: | |
panMap[iField] = iDstField | |
elif poDstLayer.CreateField( oFieldDefn ) == 0: | |
# now that we've created a field, GetLayerDefn() won't return NULL | |
if poDstFDefn is None: | |
poDstFDefn = poDstLayer.GetLayerDefn() | |
# Sanity check : if it fails, the driver is buggy | |
if poDstFDefn is not None and \ | |
poDstFDefn.GetFieldCount() != nDstFieldCount + 1: | |
print("The output driver has claimed to have added the %s field, but it did not!" % oFieldDefn.GetNameRef() ) | |
else: | |
panMap[iField] = nDstFieldCount | |
nDstFieldCount = nDstFieldCount + 1 | |
else: | |
# For an existing layer, build the map by fetching the index in the destination | |
# layer for each source field | |
if poDstFDefn is None: | |
print( "poDstFDefn == NULL.\n" ) | |
return None | |
for iField in range(nSrcFieldCount): | |
poSrcFieldDefn = poSrcFDefn.GetFieldDefn(iField) | |
iDstField = poDstFDefn.GetFieldIndex(poSrcFieldDefn.GetNameRef()) | |
if iDstField >= 0: | |
panMap[iField] = iDstField | |
iSrcZField = -1 | |
if pszZField is not None: | |
iSrcZField = poSrcFDefn.GetFieldIndex(pszZField) | |
psInfo = TargetLayerInfo() | |
psInfo.poDstLayer = poDstLayer | |
psInfo.poCT = poCT | |
#psInfo.papszTransformOptions = papszTransformOptions | |
psInfo.panMap = panMap | |
psInfo.iSrcZField = iSrcZField | |
psInfo.setFields = set() | |
psInfo.mapSrcNameToDstName = {} | |
for iField in range(poDstLayer.GetLayerDefn().GetFieldCount()): | |
name = poDstLayer.GetLayerDefn().GetFieldDefn(iField).GetNameRef() | |
psInfo.setFields.add( name ) | |
psInfo.mapSrcNameToDstName[name] = name | |
return psInfo | |
#********************************************************************** | |
# TranslateLayer() | |
#********************************************************************** | |
def TranslateLayer( psInfo, poSrcDS, poSrcLayer, poDstDS, \ | |
poOutputSRS, bNullifyOutputSRS, \ | |
eGType, bPromoteToMulti, nCoordDim, eGeomOp, dfGeomOpParam, \ | |
nCountLayerFeatures, \ | |
poClipSrc, poClipDst, bExplodeCollections, nSrcFileSize, \ | |
pnReadFeatureCount, pfnProgress, pProgressArg) : | |
bForceToPolygon = False | |
bForceToMultiPolygon = False | |
bForceToMultiLineString = False | |
poDstLayer = psInfo.poDstLayer | |
#papszTransformOptions = psInfo.papszTransformOptions | |
poCT = psInfo.poCT | |
panMap = psInfo.panMap | |
iSrcZField = psInfo.iSrcZField | |
is_shape_file = (poDstDS.GetDriver().GetName() == 'ESRI Shapefile') | |
if poOutputSRS is None and not bNullifyOutputSRS: | |
poOutputSRS = poSrcLayer.GetSpatialRef() | |
if wkbFlatten(eGType) == ogr.wkbPolygon: | |
bForceToPolygon = True | |
elif wkbFlatten(eGType) == ogr.wkbMultiPolygon: | |
bForceToMultiPolygon = True | |
elif wkbFlatten(eGType) == ogr.wkbMultiLineString: | |
bForceToMultiLineString = True | |
# -------------------------------------------------------------------- | |
# Transfer features. | |
# -------------------------------------------------------------------- | |
nFeaturesInTransaction = 0 | |
nCount = 0 | |
if nGroupTransactions > 0: | |
poDstLayer.StartTransaction() | |
while True: | |
poDstFeature = None | |
if nFIDToFetch != ogr.NullFID: | |
#// Only fetch feature on first pass. | |
if nFeaturesInTransaction == 0: | |
poFeature = poSrcLayer.GetFeature(nFIDToFetch) | |
else: | |
poFeature = None | |
else: | |
poFeature = poSrcLayer.GetNextFeature() | |
if poFeature is None: | |
break | |
nParts = 0 | |
nIters = 1 | |
if bExplodeCollections: | |
poSrcGeometry = poFeature.GetGeometryRef() | |
if poSrcGeometry is not None: | |
eSrcType = wkbFlatten(poSrcGeometry.GetGeometryType()) | |
if eSrcType == ogr.wkbMultiPoint or \ | |
eSrcType == ogr.wkbMultiLineString or \ | |
eSrcType == ogr.wkbMultiPolygon or \ | |
eSrcType == ogr.wkbGeometryCollection: | |
nParts = poSrcGeometry.GetGeometryCount() | |
nIters = nParts | |
if nIters == 0: | |
nIters = 1 | |
for iPart in range(nIters): | |
nFeaturesInTransaction = nFeaturesInTransaction + 1 | |
if nFeaturesInTransaction == nGroupTransactions: | |
poDstLayer.CommitTransaction() | |
poDstLayer.StartTransaction() | |
nFeaturesInTransaction = 0 | |
idx = poFeature.GetFieldIndex('other_tags') | |
tags = [] | |
if idx >= 0 and poFeature.IsFieldSet(idx): | |
other_tags = poFeature.GetFieldAsString(idx) | |
i = 0 | |
length = len(other_tags) | |
while True: | |
assert(other_tags[i] == '"') | |
i+= 1 | |
key = '' | |
while True: | |
if other_tags[i] == '\\': | |
key += other_tags[i+1] | |
i += 2 | |
elif other_tags[i] == '"': | |
break | |
else: | |
key += other_tags[i] | |
i += 1 | |
assert(other_tags[i] == '"') | |
i += 1 | |
assert(other_tags[i] == '=') | |
i += 1 | |
assert(other_tags[i] == '>') | |
i += 1 | |
assert(other_tags[i] == '"') | |
i+= 1 | |
value = '' | |
while True: | |
if other_tags[i] == '\\': | |
value += other_tags[i+1] | |
i += 2 | |
elif other_tags[i] == '"': | |
break | |
else: | |
value += other_tags[i] | |
i += 1 | |
assert(other_tags[i] == '"') | |
i += 1 | |
#print(key, value) | |
if not key.startswith('traffic:hourly:') and not key.startswith('seamark:light:'): | |
tags += [ (key, value) ] | |
if i == length: | |
break | |
assert(other_tags[i] == ',') | |
i += 1 | |
for (tag_name,v) in tags: | |
if tag_name not in psInfo.setFields: | |
fld_defn = ogr.FieldDefn( tag_name, ogr.OFTString) | |
if is_shape_file: | |
fld_defn.SetWidth( len(v) ) | |
poDstLayer.CreateField( fld_defn ) | |
psInfo.setFields.add( tag_name ) | |
psInfo.mapSrcNameToDstName[tag_name] = poDstLayer.GetLayerDefn().GetFieldDefn( poDstLayer.GetLayerDefn().GetFieldCount() - 1 ).GetNameRef() | |
gdal.ErrorReset() | |
poDstFeature = ogr.Feature( poDstLayer.GetLayerDefn() ) | |
if poDstFeature.SetFromWithMap( poFeature, 1, panMap ) != 0: | |
if nGroupTransactions > 0: | |
poDstLayer.CommitTransaction() | |
print("Unable to translate feature %d from layer %s" % (poFeature.GetFID() , poSrcLayer.GetName() )) | |
return False | |
for (tag_name,v) in tags: | |
poDstFeature.SetField( psInfo.mapSrcNameToDstName[tag_name], v ) | |
if bPreserveFID: | |
poDstFeature.SetFID( poFeature.GetFID() ) | |
poDstGeometry = poDstFeature.GetGeometryRef() | |
if poDstGeometry is not None: | |
if nParts > 0: | |
# For -explodecollections, extract the iPart(th) of the geometry | |
poPart = poDstGeometry.GetGeometryRef(iPart).Clone() | |
poDstFeature.SetGeometryDirectly(poPart) | |
poDstGeometry = poPart | |
if iSrcZField != -1: | |
SetZ(poDstGeometry, poFeature.GetFieldAsDouble(iSrcZField)) | |
# This will correct the coordinate dimension to 3 | |
poDupGeometry = poDstGeometry.Clone() | |
poDstFeature.SetGeometryDirectly(poDupGeometry) | |
poDstGeometry = poDupGeometry | |
if nCoordDim == 2 or nCoordDim == 3: | |
poDstGeometry.SetCoordinateDimension( nCoordDim ) | |
if eGeomOp == GeomOperation.SEGMENTIZE: | |
pass | |
#if (poDstFeature.GetGeometryRef() is not None and dfGeomOpParam > 0) | |
# poDstFeature.GetGeometryRef().segmentize(dfGeomOpParam); | |
elif eGeomOp == GeomOperation.SIMPLIFY_PRESERVE_TOPOLOGY and dfGeomOpParam > 0: | |
poNewGeom = poDstGeometry.SimplifyPreserveTopology(dfGeomOpParam) | |
if poNewGeom is not None: | |
poDstFeature.SetGeometryDirectly(poNewGeom) | |
poDstGeometry = poNewGeom | |
if poClipSrc is not None: | |
poClipped = poDstGeometry.Intersection(poClipSrc) | |
if poClipped is None or poClipped.IsEmpty(): | |
# Report progress | |
nCount = nCount +1 | |
if pfnProgress is not None: | |
pfnProgress(nCount * 1.0 / nCountLayerFeatures, "", pProgressArg) | |
continue | |
poDstFeature.SetGeometryDirectly(poClipped) | |
poDstGeometry = poClipped | |
if poCT is not None: | |
eErr = poDstGeometry.Transform( poCT ) | |
if eErr != 0: | |
if nGroupTransactions > 0: | |
poDstLayer.CommitTransaction() | |
print("Failed to reproject feature %d (geometry probably out of source or destination SRS)." % poFeature.GetFID()) | |
if not bSkipFailures: | |
return False | |
elif poOutputSRS is not None: | |
poDstGeometry.AssignSpatialReference(poOutputSRS) | |
if poClipDst is not None: | |
poClipped = poDstGeometry.Intersection(poClipDst) | |
if poClipped is None or poClipped.IsEmpty(): | |
continue | |
poDstFeature.SetGeometryDirectly(poClipped) | |
poDstGeometry = poClipped | |
if bForceToPolygon: | |
poDstFeature.SetGeometryDirectly(ogr.ForceToPolygon(poDstGeometry)) | |
elif bForceToMultiPolygon or \ | |
(bPromoteToMulti and wkbFlatten(poDstGeometry.GetGeometryType()) == ogr.wkbPolygon): | |
poDstFeature.SetGeometryDirectly(ogr.ForceToMultiPolygon(poDstGeometry)) | |
elif bForceToMultiLineString or \ | |
(bPromoteToMulti and wkbFlatten(poDstGeometry.GetGeometryType()) == ogr.wkbLineString): | |
poDstFeature.SetGeometryDirectly(ogr.ForceToMultiLineString(poDstGeometry)) | |
gdal.ErrorReset() | |
if poDstLayer.CreateFeature( poDstFeature ) != 0 and not bSkipFailures: | |
if nGroupTransactions > 0: | |
poDstLayer.RollbackTransaction() | |
return False | |
# Report progress | |
nCount = nCount + 1 | |
if pfnProgress is not None: | |
if nSrcFileSize != 0: | |
if (nCount % 1000) == 0: | |
poFCLayer = poSrcDS.ExecuteSQL("GetBytesRead()", None, None) | |
if poFCLayer is not None: | |
poFeat = poFCLayer.GetNextFeature() | |
if poFeat is not None: | |
pszReadSize = poFeat.GetFieldAsString(0) | |
nReadSize = int(pszReadSize) | |
pfnProgress(nReadSize * 1.0 / nSrcFileSize, "", pProgressArg) | |
poSrcDS.ReleaseResultSet(poFCLayer) | |
else: | |
pfnProgress(nCount * 1.0 / nCountLayerFeatures, "", pProgressArg) | |
if pnReadFeatureCount is not None: | |
pnReadFeatureCount[0] = nCount | |
if nGroupTransactions > 0: | |
poDstLayer.CommitTransaction() | |
return True | |
if __name__ == '__main__': | |
version_num = int(gdal.VersionInfo('VERSION_NUM')) | |
if version_num < 1800: # because of ogr.GetFieldTypeName | |
print('ERROR: Python bindings of GDAL 1.8.0 or later required') | |
sys.exit(1) | |
if not main(sys.argv): | |
sys.exit(1) | |
else: | |
sys.exit(0) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment