#!/usr/bin/env python3 # -*- coding: utf-8 -*- ############################################################################### # $Id$ # # Project: OGR Python samples # Purpose: Create OGR VRT from source datasource # Author: Frank Warmerdam, warmerdam@pobox.com # ############################################################################### # Copyright (c) 2009, Frank Warmerdam # Copyright (c) 2009-2014, Even Rouault # # 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.path import sys from osgeo import gdal, ogr ############################################################################# def GeomType2Name(typ): flat_type = ogr.GT_Flatten(typ) dic = { ogr.wkbUnknown: ("wkbUnknown", "25D"), ogr.wkbPoint: ("wkbPoint", "25D"), ogr.wkbLineString: ("wkbLineString", "25D"), ogr.wkbPolygon: ("wkbPolygon", "25D"), ogr.wkbMultiPoint: ("wkbMultiPoint", "25D"), ogr.wkbMultiLineString: ("wkbMultiLineString", "25D"), ogr.wkbMultiPolygon: ("wkbMultiPolygon", "25D"), ogr.wkbGeometryCollection: ("wkbGeometryCollection", "25D"), ogr.wkbNone: ("wkbNone", ""), ogr.wkbLinearRing: ("wkbLinearRing", ""), ogr.wkbCircularString: ("wkbCircularString", "Z"), ogr.wkbCompoundCurve: ("wkbCompoundCurve", "Z"), ogr.wkbCurvePolygon: ("wkbCurvePolygon", "Z"), ogr.wkbMultiCurve: ("wkbMultiCurve", "Z"), ogr.wkbMultiSurface: ("wkbMultiSurface", "Z"), ogr.wkbCurve: ("wkbCurve", "Z"), ogr.wkbSurface: ("wkbSurface", "Z"), ogr.wkbPolyhedralSurface: ("wkbPolyhedralSurface", "Z"), ogr.wkbTIN: ("wkbTIN", "Z"), ogr.wkbTriangle: ("wkbTriangle", "Z"), } ret = dic[flat_type][0] if flat_type != typ: if ogr.GT_HasM(typ): if ogr.GT_HasZ(typ): ret += "ZM" else: ret += "M" else: ret += dic[flat_type][1] return ret ############################################################################# def Esc(x): return gdal.EscapeString(x, gdal.CPLES_XML) ############################################################################# def Usage(): print("Usage: ogr2vrt.py [-relative] [-schema] [-feature_count] [-extent]") print(" in_datasource out_vrtfile [layers]") print("") return 2 def main(argv=sys.argv): infile = None outfile = None layer_list = [] relative = "0" schema = 0 feature_count = 0 extent = 0 openoptions = [] argv = gdal.GeneralCmdLineProcessor(argv) if argv is None: return 0 i = 1 while i < len(argv): arg = argv[i] if arg == "-relative": relative = "1" elif arg == "-schema": schema = 1 elif arg == "-feature_count": feature_count = 1 elif arg == "-extent": extent = 1 elif arg == "-oo": i += 1 openoptions.append(argv[i]) elif arg[0] == "-": return Usage() elif infile is None: infile = arg elif outfile is None: outfile = arg else: layer_list.append(arg) i = i + 1 if outfile is None: return Usage() if schema and feature_count: sys.stderr.write("Ignoring -feature_count when used with -schema.\n") feature_count = 0 if schema and extent: sys.stderr.write("Ignoring -extent when used with -schema.\n") extent = 0 ############################################################################# # Open the datasource to read. src_ds = gdal.OpenEx(infile, gdal.OF_VECTOR, open_options=openoptions) if schema: infile = "@dummy@" if not layer_list: for lyr_idx in range(src_ds.GetLayerCount()): layer_list.append(src_ds.GetLayer(lyr_idx).GetLayerDefn().GetName()) ############################################################################# # Start the VRT file. vrt = "\n" ############################################################################# # Metadata mdd_list = src_ds.GetMetadataDomainList() if mdd_list is not None: for domain in mdd_list: if domain == "": vrt += " \n" elif len(domain) > 4 and domain[0:4] == "xml:": vrt += ' \n' % Esc(domain) else: vrt += ' \n' % Esc(domain) if len(domain) > 4 and domain[0:4] == "xml:": vrt += src_ds.GetMetadata_List(domain)[0] else: md = src_ds.GetMetadata(domain) for key in md: vrt += ' %s\n' % (Esc(key), Esc(md[key])) vrt += " \n" ############################################################################# # Process each source layer. for name in layer_list: layer = src_ds.GetLayerByName(name) layerdef = layer.GetLayerDefn() vrt += ' \n' % Esc(name) mdd_list = layer.GetMetadataDomainList() if mdd_list is not None: for domain in mdd_list: if domain == "": vrt += " \n" elif len(domain) > 4 and domain[0:4] == "xml:": vrt += ' \n' % Esc(domain) else: vrt += ' \n' % Esc(domain) if len(domain) > 4 and domain[0:4] == "xml:": vrt += layer.GetMetadata_List(domain)[0] else: md = layer.GetMetadata(domain) for key in md: vrt += ' %s\n' % ( Esc(key), Esc(md[key]), ) vrt += " \n" if ( not os.path.isabs(outfile) and not os.path.isabs(infile) and os.path.dirname(outfile) == "" and os.path.dirname(infile) == "" ): relative = 1 vrt += ( ' %s\n' % (relative, not schema, Esc(infile)) ) if openoptions: vrt += " \n" for option in openoptions: (key, value) = option.split("=") vrt += ' %s\n' % (Esc(key), Esc(value)) vrt += " \n" if schema: vrt += " @dummy@\n" else: vrt += " %s\n" % Esc(name) # Historic format for mono-geometry layers if layerdef.GetGeomFieldCount() == 0: vrt += " wkbNone\n" elif ( layerdef.GetGeomFieldCount() == 1 and layerdef.GetGeomFieldDefn(0).IsNullable() ): vrt += " %s\n" % GeomType2Name( layerdef.GetGeomType() ) srs = layer.GetSpatialRef() if srs is not None: vrt += " %s\n" % (Esc(srs.ExportToWkt())) if extent: (xmin, xmax, ymin, ymax) = layer.GetExtent() vrt += " %.15g\n" % xmin vrt += " %.15g\n" % ymin vrt += " %.15g\n" % xmax vrt += " %.15g\n" % ymax # New format for multi-geometry field support else: for fld_index in range(layerdef.GetGeomFieldCount()): src_fd = layerdef.GetGeomFieldDefn(fld_index) vrt += ' 0: vrt += ' width="%d"' % src_fd.GetWidth() if src_fd.GetPrecision() > 0: vrt += ' precision="%d"' % src_fd.GetPrecision() if src_fd.IsNullable() == 0: vrt += ' nullable="false"' if src_fd.GetAlternativeName(): vrt += ' alternativeName="%s"' % Esc(src_fd.GetAlternativeName()) if src_fd.GetComment(): vrt += ' comment="%s"' % Esc(src_fd.GetComment()) try: if src_fd.IsUnique(): vrt += ' unique="true"' except AttributeError: # if run with GDAL < 3.2 pass vrt += "/>\n" if feature_count: vrt += " %d\n" % layer.GetFeatureCount() vrt += " \n" vrt += "\n" ############################################################################# # Write vrt open(outfile, "w").write(vrt) return 0 if __name__ == "__main__": sys.exit(main(sys.argv))