EVOLUTION-MANAGER

Edit File: densify.py

#!/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()