source: subversion/applications/utils/osm-extract/polygons/ogr2poly.py @ 29411

Last change on this file since 29411 was 27188, checked in by joshdoe, 8 years ago

add ogr2poly.py which converts OGR supported files to POLY files

File size: 8.2 KB
Line 
1#!/usr/bin/env python
2
3# This converts OGR supported files (Shapefile, GPX, etc.) to the polygon
4# filter file format [1] supported by Osmosis and other tools. If there is
5# more than one feature, it will create one POLY file for each feature,
6# either using an incrementing filename or based on a field value. It also
7# includes buffering and simplifying. This allows point or line features
8# to be used when creating POLY files, but in this case buffering must
9# be used.
10#
11# Requires GDAL/OGR compiled with GEOS
12#
13# [1] http://wiki.openstreetmap.org/wiki/Osmosis/Polygon_Filter_File_Format
14#
15# written by Josh Doe <josh@joshdoe.com> and licensed under the LGPL
16#
17# This program is free software: you can redistribute it and/or modify
18# it under the terms of the GNU Lesser General Public License as published by
19# the Free Software Foundation, either version 3 of the License, or
20# (at your option) any later version.
21
22# This program is distributed in the hope that it will be useful,
23# but WITHOUT ANY WARRANTY; without even the implied warranty of
24# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
25# GNU Lesser General Public License for more details.
26
27# You should have received a copy of the GNU Lesser General Public License
28# along with this program.  If not, see <http://www.gnu.org/licenses/>.
29
30from optparse import OptionParser
31import logging
32import os
33import sys
34
35from osgeo import ogr
36from osgeo import osr
37
38# TODO:
39#  check if file exists, make sure field is unique (increment)
40#  likely doesn't handle areas spanning the antimeridian (+/-180 deg lon)
41#  sometimes an empty poly is created, so pay attention to warnings
42#      this can usually be fixed by decreasing the simplify distance
43
44
45def createPolys(inOgr, options):
46    logging.info("Opening datasource '%s'" % inOgr)
47    ds = ogr.Open(inOgr)
48    lyr = ds.GetLayer(options.layer)
49
50    # create SRS transformations
51    mercSRS = osr.SpatialReference()
52    mercSRS.ImportFromEPSG(3857)  # TODO: make this an option
53    wgsSRS = osr.SpatialReference()
54    wgsSRS.ImportFromEPSG(4326)
55    nativeSRS2bufferSRS = osr.CoordinateTransformation(lyr.GetSpatialRef(),
56                                                       mercSRS)
57    bufferSRS2wgsSRS = osr.CoordinateTransformation(mercSRS, wgsSRS)
58    nativeSRS2wgsSRS = osr.CoordinateTransformation(lyr.GetSpatialRef(),
59                                                    wgsSRS)
60
61    # if no field name is provided, use incrementing number
62    # (padded with just enough zeros)
63    inc = 0
64    incFmt = '%0' + str(len(str(lyr.GetFeatureCount() - 1))) + 'd'
65
66    logging.info('Found %d features, will create one POLY file for each one'
67           % lyr.GetFeatureCount())
68
69    # create POLYs
70    for feat in lyr:
71        if options.fieldName != None:
72            fieldVal = feat.GetFieldAsString(options.fieldName)
73            if fieldVal is None:
74                return False
75            polyName = options.outPrefix + fieldVal.replace(' ', '_')
76        else:
77            polyName = options.outPrefix + incFmt % inc
78            inc += 1
79
80        logging.info('Creating ' + polyName + '.poly')
81        f = open(polyName + '.poly', 'wt')
82        print >>f, polyName
83
84        # this will be a polygon, TODO: handle linestrings (must be buffered)
85        geom = feat.GetGeometryRef()
86        geomType = geom.GetGeometryType()
87
88        subGeom = []
89
90        nonAreaTypes = [ogr.wkbPoint, ogr.wkbLineString, ogr.wkbMultiPoint,
91                        ogr.wkbMultiLineString]
92        if geomType in nonAreaTypes and options.bufferDistance == 0:
93            logging.warn("Ignoring non-area type. " +
94                         "To include you must set a buffer distance.")
95            continue
96        if geomType in [ogr.wkbUnknown, ogr.wkbNone]:
97            logging.warn("Ignoring unknown geometry type.")
98            continue
99
100        # transform to WGS84, buffering/simplifying along the way
101        if options.bufferDistance > 0 or options.simplifyDistance > 0:
102            geom.Transform(nativeSRS2bufferSRS)
103
104            if options.bufferDistance > 0:
105                geom = geom.Buffer(float(options.bufferDistance))
106            if options.simplifyDistance > 0:
107                geom = geom.Simplify(float(options.simplifyDistance))
108
109            geom.Transform(bufferSRS2wgsSRS)
110        else:
111            geom.Transform(nativeSRS2wgsSRS)
112
113        # handle multi-polygons
114        subgeom = []
115        geomtype = geom.GetGeometryType()
116        if geomtype == ogr.wkbPolygon:
117            subgeom = [geom]
118        elif geomtype == ogr.wkbMultiPolygon:
119            for k in range(geom.GetGeometryCount()):
120                subgeom.append(geom.GetGeometryRef(k))
121
122        logging.debug("# of polygons: " + str(len(subgeom)))
123        for g in subgeom:
124            # loop over all rings in the polygon
125            logging.debug('# of rings: ' + str(g.GetGeometryCount()))
126            for i in range(0, g.GetGeometryCount()):
127                if i == 0:
128                    # outer ring
129                    print >>f, i + 1
130                else:
131                    # inner ring
132                    print >>f, '!%d' % (i + 1)
133                ring = g.GetGeometryRef(i)
134
135                if ring.GetPointCount() > 0:
136                    logging.debug('# of points: ' + str(ring.GetPointCount()))
137                else:
138                    logging.warn('Ring with no points')
139
140                # output all points in the ring
141                for j in range(0, ring.GetPointCount()):
142                    (x, y, z) = ring.GetPoint(j)
143                    print >>f, '   %.6E   %.6E' % (x, y)
144                print >>f, 'END'
145        print >>f, 'END'
146        f.close()
147    return True
148
149if __name__ == '__main__':
150    # Setup program usage
151    usage = "Usage: %prog [options] src_datasource_name [layer]"
152    parser = OptionParser(usage=usage)
153    parser.add_option("-p", "--prefix", dest="outPrefix",
154                      help="Text to prepend to POLY filenames.")
155    parser.add_option("-b", "--buffer-distance", dest="bufferDistance",
156                      type="float",
157                      help="Set buffer distance in meters (default: 0).")
158    parser.add_option("-s", "--simplify-distance", dest="simplifyDistance",
159                      type="float",
160                      help="Set simplify tolerance in meters (default: 0).")
161    parser.add_option("-f", "--field-name", dest="fieldName",
162                      help="Field name to use to name files.")
163    parser.add_option("-v", "--verbose", dest="verbose", action="store_true",
164                      help="Print detailed status messages.")
165
166    parser.set_defaults(bufferDistance=0, fieldName=None, outPrefix=None,
167        simplifyDistance=0, layer=0, verbose=False)
168
169    # Parse and process arguments
170    (options, args) = parser.parse_args()
171
172    if options.verbose:
173        logging.basicConfig(format='%(asctime)s:%(levelname)s: %(message)s',
174                            level=logging.DEBUG)
175    else:
176        logging.basicConfig(format='%(levelname)s: %(message)s',
177                            level=logging.WARNING)
178
179    if len(args) < 1:
180        parser.print_help()
181        parser.error("You must specify an OGR source")
182        sys.exit(1)
183    elif len(args) > 2:
184        parser.error("You have specified too many arguments")
185
186    # note that this may be a file (e.g. .shp) or a database connection string
187    src_datasource = args[0]
188    if len(args) == 2:
189        options.layer = args[1]
190
191    # check options
192    if options.outPrefix == None:
193        if os.path.exists(src_datasource):
194            # put in current dir, TODO: allow user to specify output dir?
195            (options.outPrefix, ext) = os.path.splitext(
196                    os.path.basename(src_datasource))
197            options.outPrefix += '_'
198        else:
199            # file doesn't exist, so possibly a DB connection string
200            options.outPrefix = 'poly_'
201    if options.bufferDistance < 0:
202        parser.error("Buffer distance must be greater than zero.")
203    if options.simplifyDistance < 0:
204        parser.error("Simplify tolerance must be greater than zero.")
205    if options.simplifyDistance > options.bufferDistance:
206        logging.warn("Simplify distance greater than buffer distance")
207
208    if createPolys(src_datasource, options):
209        logging.info('Finished!')
210        sys.exit(0)
211    else:
212        logging.info('Failed!')
213        sys.exit(1)
Note: See TracBrowser for help on using the repository browser.