File: densify.py

package info (click to toggle)
gdal 1.10.1%2Bdfsg-8
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 84,320 kB
  • ctags: 74,726
  • sloc: cpp: 677,199; ansic: 162,820; python: 13,816; cs: 11,163; sh: 10,446; java: 5,279; perl: 4,429; php: 2,971; xml: 1,500; yacc: 934; makefile: 494; sql: 112
file content (381 lines) | stat: -rwxr-xr-x 13,884 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
#!/usr/bin/env python

try:
    from osgeo import osr
    from osgeo import ogr
    ogr.UseExceptions()
except ImportError:
    import osr
    import ogr

import sys
import os
import math




class Translator(object):

    def construct_parser(self):
        from optparse import OptionParser, OptionGroup
        usage = "usage: %prog [options] arg"
        parser = OptionParser(usage)
        g = OptionGroup(parser, "Base options", "Basic Translation Options")
        g.add_option("-i", "--input", dest="input",
                          help="OGR input data source", metavar="INPUT")
        g.add_option("-o", "--output", dest='output',
                          help="OGR output data source", metavar="OUTPUT")
        g.add_option("-n", "--nooverwrite",
                          action="store_false", dest="overwrite", 
                          help="Do not overwrite the existing output file")

        g.add_option("-f", "--driver", dest='driver',
                          help="OGR output driver.  Defaults to \"ESRI Shapefile\"", metavar="DRIVER")

        g.add_option('-w', "--where", dest="where",
                          help="""SQL attribute filter -- enclose in quotes "PNAME='Cedar River'" """,)
        g.add_option('-s', "--spat", dest="spat",
                          help="""Spatial query extents -- minx miny maxx maxy""",
                          type=float, nargs=4)
        g.add_option('-l', "--layer", dest="layer",
                          help="""The name of the input layer to translate, if not given, the first layer on the data source is used""",)                          

        g.add_option('-k', "--select", dest="fields",
                          help="""Comma separated list of fields to include -- field1,field2,field3,...,fieldn""",)                          
        g.add_option('-t', '--target-srs', dest="t_srs",
                          help="""Target SRS -- the spatial reference system to project the data to""")
        g.add_option('-a', '--assign-srs', dest="a_srs",
                          help="""Assign SRS -- the spatial reference to assign to the input data""")
        g.add_option("-q", "--quiet",
                          action="store_false", dest="verbose", default=False,
                          help="don't print status messages to stdout")
        parser.add_option_group(g)

        if self.opts:
            g = OptionGroup(parser, "Special Options", "Special options")
            for o in self.opts:
                g.add_option(o)
            parser.add_option_group(g)

        parser.set_defaults(verbose=True, driver="ESRI Shapefile", overwrite=True)

        self.parser = parser
        
    def __init__(self, arguments, options=None):
        self.input = None
        self.output = None
        
        self.opts = options
        self.construct_parser()
        self.options, self.args = self.parser.parse_args(args=arguments)

    def process(self):
        self.open()
        self.make_fields()
        self.translate()
        
    def open(self):
        if self.options.input:
            self.in_ds = ogr.Open(self.options.input)
        else:
            raise Exception("No input layer was specified")
        if self.options.layer:
            self.input = self.in_ds.GetLayerByName(self.options.layer)
            if not self.input:
                raise Exception("The layer '%s' was not found" % self.options.layer)
        else:
            self.input = self.in_ds.GetLayer(0)

        if self.options.a_srs:
            self.in_srs = osr.SpatialReference()
            self.in_srs.SetFromUserInput(self.options.a_srs)
        else:
            self.in_srs = self.input.GetSpatialRef()

        if self.options.spat:
            self.input.SetSpatialFilterRect(*self.options.spat)
            
        self.out_drv = ogr.GetDriverByName(self.options.driver)
        
        if self.options.where:
            self.input.SetAttributeFilter(self.options.where)
            
        if not self.out_drv:
            raise Exception("The '%s' driver was not found, did you misspell it or is it not available in this GDAL build?", self.options.driver)
        if not self.out_drv.TestCapability( 'CreateDataSource' ):
            raise Exception("The '%s' driver does not support creating layers, you will have to choose another output driver", self.options.driver)
        if not self.options.output:
            raise Exception("No output layer was specified")
        if self.options.driver == 'ESRI Shapefile':
            path, filename = os.path.split(os.path.abspath(self.options.output))
            name, ext = os.path.splitext(filename)
            if self.options.overwrite:
            # special case the Shapefile driver, which behaves specially.
                if os.path.exists(os.path.join(path,name,) +'.shp'):
                    os.remove(os.path.join(path,name,) +'.shp')
                    os.remove(os.path.join(path,name,) +'.shx')
                    os.remove(os.path.join(path,name,) +'.dbf')
            else:
                if os.path.exists(os.path.join(path,name,)+".shp"):
                    raise Exception("The file '%s' already exists, but the overwrite option is not specified" % (os.path.join(path,name,)+".shp"))

        if self.options.overwrite:
            dsco = ('OVERWRITE=YES',)
        else:
            dsco = (),

        self.out_ds = self.out_drv.CreateDataSource( self.options.output, dsco)
        
        if self.options.t_srs:
            self.out_srs = osr.SpatialReference()
            self.out_srs.SetFromUserInput(self.options.t_srs)
        else:
            self.out_srs = None
        self.output = self.out_ds.CreateLayer(  self.options.output,
                                                geom_type = self.input.GetLayerDefn().GetGeomType(), 
                                                srs= self.out_srs)


    def make_fields(self):
        defn = self.input.GetLayerDefn()

        if self.options.fields:
            fields = self.options.fields.split(',')
            for i in range(defn.GetFieldCount()):
                fld = defn.GetFieldDefn(i)
                for f in fields:
                    if fld.GetName().upper() == f.upper():
                        self.output.CreateField(fld)
        else:
            for i in range(defn.GetFieldCount()):
                fld = defn.GetFieldDefn(i)
                self.output.CreateField(fld)

    def translate(self, geometry_callback=None, attribute_callback=None):
        f = self.input.GetNextFeature()
        trans = None
        if self.options.t_srs:
            trans = osr.CoordinateTransformation(self.in_srs, self.out_srs)
        while f:
            geom = f.GetGeometryRef().Clone()
            
            if trans:
                geom.Transform(trans)

            if geometry_callback:
                geom = geometry_callback(geom)

                
            f.SetGeometry(geom)
            d = ogr.Feature(feature_def=self.output.GetLayerDefn())
            d.SetFrom(f)
            self.output.CreateFeature(d)
            f = self.input.GetNextFeature()

            
    def __del__(self):
        if self.output:
            self.output.SyncToDisk()

def radians(degrees):
    return math.pi/180.0*degrees
def degrees(radians):
    return radians*180.0/math.pi
    
class Densify(Translator):

    def calcpoint(self,x0, x1, y0, y1, d):
        a = x1 - x0
        b = y1 - y0
        
        if a == 0:
            xn = x1
            
            if b > 0:
                yn = y0 + d
            else:
                yn = y0 - d
            return (xn, yn)
                      
        theta = degrees(math.atan(abs(b)/abs(a)))
        
        if a > 0 and b > 0:
            omega = theta
        if a < 0 and b > 0:
            omega = 180 - theta
        if a < 0 and b < 0:
            omega = 180 + theta
        if a > 0 and b < 0:
            omega = 360 - theta

        if b == 0:
            yn = y1
            if a > 0:
                xn = x0 + d
            else:
                xn = x0 - d
        else:
            xn = x0 + d*math.cos(radians(omega))
            yn = y0 + d*math.sin(radians(omega))
        
        return (xn, yn)
                    
    def distance(self, x0, x1, y0, y1):
        deltax = x0 - x1
        deltay = y0 - y1
        d2 = (deltax)**2 + (deltay)**2
        d = math.sqrt(d2)
        return d
    
    def densify(self, geometry):
        gtype = geometry.GetGeometryType()
        if  not (gtype == ogr.wkbLineString or gtype == ogr.wkbMultiLineString):
            raise Exception("The densify function only works on linestring or multilinestring geometries")
            
        g = ogr.Geometry(ogr.wkbLineString)

        # add the first point
        x0 = geometry.GetX(0)
        y0 = geometry.GetY(0)
        g.AddPoint(x0, y0)

        for i in range(1,geometry.GetPointCount()):
            threshold = self.options.distance
            x1 = geometry.GetX(i)
            y1 = geometry.GetY(i)
            if not x0 or not y0:
                raise Exception("First point is null")
            d = self.distance(x0, x1, y0, y1)

            if self.options.remainder.upper() == "UNIFORM":
                if d != 0.0:
                    threshold = float(d)/math.ceil(d/threshold)
                else:
                    # duplicate point... throw it out
                    continue
            if (d > threshold):
                if self.options.remainder.upper() == "UNIFORM":
                    segcount = int(math.ceil(d/threshold))

                    dx = (x1 - x0)/segcount
                    dy = (y1 - y0)/segcount

                    x = x0
                    y = y0
                    for p in range(1,segcount):
                        x = x + dx
                        y = y + dy
                        g.AddPoint(x, y)
                        
                elif self.options.remainder.upper() == "END":
                    segcount = int(math.floor(d/threshold))
                    xa = None
                    ya = None
                    for p in range(1,segcount):
                        if not xa:
                            xn, yn = self.calcpoint(x0,x1,y0,y1,threshold)
                            d = self.distance(x0, xn, y0, yn)
                            xa = xn
                            ya = yn
                            g.AddPoint(xa,ya)
                            continue
                        xn, yn = self.calcpoint(xa, x1, ya, y1, threshold)
                        xa = xn
                        ya = yn
                        g.AddPoint(xa,ya)
                        
                elif self.options.remainder.upper() == "BEGIN":
                    
                    # I think this might put an extra point in at the end of the 
                    # first segment
                    segcount = int(math.floor(d/threshold))
                    xa = None
                    ya = None
                    xb = x0
                    yb = y0
                    remainder = d % threshold
                    for p in range(segcount):
                        if not xa:
                            xn, yn = self.calcpoint(x0,x1,y0,y1,remainder)
 
                            d = self.distance(x0, xn, y0, yn)
                            xa = xn
                            ya = yn
                            g.AddPoint(xa,ya)
                            continue
                        xn, yn = self.calcpoint(xa, x1, ya, y1, threshold)
                        xa = xn
                        ya = yn
                        g.AddPoint(xa,ya)

            g.AddPoint(x1,y1)
            x0 = x1
            y0 = y1

                
        return g

    def process(self):
        self.open()
        self.make_fields()
        self.translate(geometry_callback = self.densify)

def GetLength(geometry):

    def get_distance(x1, y1, x2, y2):
        """Return the euclidian distance between this point and another"""
        import math
        deltax = x1 - x2
        deltay = y1 - y2
        d2 = (deltax**2) + (deltay**2)
        distance = math.sqrt(d2)
        return distance

    def cumulate(single):
        cumulative = 0.0
        g = single
        pt_count = g.GetPointCount()
        x1, y1 = g.GetX(0), g.GetY(0)
        for pi in range(1,pt_count):
             x2, y2 = g.GetX(pi), g.GetY(pi)
             length = get_distance(x1, y1, x2, y2)
             cumulative = cumulative + length
             x1 = x2
             y1 = y2
        return cumulative

    cumulative = 0.0
    geom_count = geometry.GetGeometryCount()
    
    if geom_count:
        for gi in range(geom_count):
            g = geometry.GetGeometryRef(gi)
            cumulative = cumulative + cumulate(g)
    else:
        cumulative = cumulate(geometry)
    return cumulative
    
def main():
    import optparse

    options = []
    o = optparse.make_option("-r", "--remainder", dest="remainder",
                         type="choice",default='end', 
                          help="""what to do with the remainder -- place it at the beginning, 
place it at the end, or evenly distribute it across the segment""",
                          choices=['end','begin','uniform'])
    options.append(o)
    o = optparse.make_option("-d", "--distance", dest='distance', type="float",
                          help="""Threshold distance for point placement.  If the 
'uniform' remainder is used, points will be evenly placed 
along the segment in a fashion that makes sure they are 
no further apart than the threshold.  If 'beg' or 'end' 
is chosen, the threshold distance will be used as an absolute value.""", 
                          metavar="DISTANCE")
    options.append(o)
    d = Densify(sys.argv[1:], options=options)
    d.process()

if __name__=='__main__':
    main()