I've tried to adapt a python ogr script which displaces points too near to each other (originally from http://wiki.gis-lab.info/w/Displacement_of_points_with_same_coordinates,_ShiftPoints/OGR for python2) to Python3. I've also changed the GDAL driver from Shapefile to GeoJSON:
#!/usr/bin/env python
import sys
import math
from itertools import cycle
from optparse import OptionParser
try:
from osgeo import ogr
except ImportError:
import ogr
driverName = "GeoJSON"
class progressBar( object ):
'''Class to display progress bar
'''
def __init__( self, maximum, barLength ):
'''Init progressbar instance.
@param maximum maximum progress value
@param barLength length of the bar in characters
'''
self.maxValue = maximum
self.barLength = barLength
self.spin = cycle(r'-\|/').__next__
self.lastLength = 0
self.tmpl = '%-' + str( barLength ) + 's ] %c %5.1f%%'
sys.stdout.write( '[ ' )
sys.stdout.flush()
def update( self, value ):
'''Update progressbar.
@param value Input new progress value
'''
# Remove last state.
sys.stdout.write( '\b' * self.lastLength )
percent = value * 100.0 / self.maxValue
# Generate new state
width = int( percent / 100.0 * self.barLength )
output = self.tmpl % ( '-' * width, self.spin(), percent )
# Show the new state and store its length.
sys.stdout.write( output )
sys.stdout.flush()
self.lastLength = len( output )
def __log( msg, abort = True, exitCode = 1 ):
''' display log message and abort execution
@param msg message to display
@param abort abort program execution or no. Default: True
@param exitCode exit code passed to sys.exit. Default: 1
'''
print(msg)
if abort:
sys.exit( exitCode )
def doDisplacement( src, dst, radius, rotate, logFile = None, logField = None ):
''' read source file and write displacement result to output file
@param src source shapefile
@param dst destination shapefile
@param radius displacement radius
@param rotate flag indicates rotate points or no
'''
# try to open source shapefile
inShape = ogr.Open( src, False )
if inShape is None:
__log( "Unable to open source shapefile " + src )
print(inShape)
inLayer = inShape.GetLayer( 0 )
print(inLayer)
print(inLayer.GetGeomType())
if inLayer.GetGeomType() not in [ ogr.wkbPoint, ogr.wkbPoint25D ]:
__log( "Input layer should be point layer." )
print(inLayer.GetFeatureCount())
if inLayer.GetFeatureCount() == 0:
__log( "There are no points in given shapefile." )
drv = ogr.GetDriverByName( driverName )
if drv is None:
__log( "%s driver not available." % driverName )
crs = inLayer.GetSpatialRef()
# initialize output shapefile
outShape = drv.CreateDataSource( dst )
if outShape is None:
__log( "Creation of output file %s failed." % dst )
outLayer = outShape.CreateLayer( "point", crs, ogr.wkbPoint )
if outLayer is None:
__log( "Layer creation failed." )
# create fields in output layer
logIndex = 0
featDef = inLayer.GetLayerDefn()
for i in range( featDef.GetFieldCount() ):
fieldDef = featDef.GetFieldDefn( i )
if logField and logField.lower() == fieldDef.GetNameRHef().lower():
logIndex = i
if outLayer.CreateField ( fieldDef ) != 0:
__log( "Creating %s field failed." % fieldDef.GetNameRef() )
__log( "Search for overlapped features", abort = False )
cn = 0
pb = progressBar( inLayer.GetFeatureCount(), 65 )
# find overlapped points
displacement = dict()
inLayer.ResetReading()
inFeat = inLayer.GetNextFeature()
while inFeat is not None:
wkt = inFeat.GetGeometryRef().ExportToWkt()
print(wkt)
if wkt not in displacement:
lst = [ inFeat.GetFID() ]
displacement[ wkt ] = lst
else:
lst = displacement[ wkt ]
j = [ inFeat.GetFID() ]
j.extend( lst )
displacement[ wkt ] = j
cn += 1
pb.update( cn )
inFeat = inLayer.GetNextFeature()
__log( "\nDisplace overlapped features", abort = False )
cn = 0
pb = progressBar( len( displacement ), 65 )
lineIds = []
lineValues = []
lines = []
pc = 0
# move overlapped features and write them to output file
for k, v in displacement.items():
featNum = len( v )
if featNum == 1:
f = inLayer.GetFeature( v[ 0 ] )
if outLayer.CreateFeature( f ) != 0:
__log( "Failed to create feature in output shapefile." )
else:
pc += featNum
fullPerimeter = 2 * math.pi
angleStep = fullPerimeter / featNum
if featNum == 2 and rotate:
currentAngle = math.pi / 2
else:
currentAngle = 0
for i in v:
sinusCurrentAngle = math.sin( currentAngle )
cosinusCurrentAngle = math.cos( currentAngle )
dx = radius * sinusCurrentAngle
dy = radius * cosinusCurrentAngle
ft = inLayer.GetFeature( i )
lineIds.append( str( ft.GetFID() ) )
lineValues.append( ft.GetFieldAsString( logIndex ) )
geom = ft.GetGeometryRef()
x = geom.GetX()
y = geom.GetY()
pt = ogr.Geometry( ogr.wkbPoint )
pt.SetPoint_2D( 0, x + dx, y + dy )
f = ogr.Feature( outLayer.GetLayerDefn() )
f.SetFrom( ft )
f.SetGeometry( pt )
if outLayer.CreateFeature( f ) != 0:
__log( "Failed to create feature in output shapefile." )
f.Destroy()
currentAngle += angleStep
# store logged info
lines.append( ";".join( lineIds ) + "(" + ";".join( lineValues ) + ")\n" )
lineIds = []
lineValues = []
cn += 1
pb.update( cn )
# cleanup
print("\nPoints displaced: %d\n" % pc)
inShape = None
outShape = None
# write log
if logFile:
fl = open( logFile, "w" )
fl.writelines( lines )
fl.close()
if __name__ == '__main__':
parser = OptionParser( usage = "displacement.py [OPTIONS] INPUT OUTPUT",
version = "%prog v0.1",
description = "Shift overlapped points so all of them becomes visible",
epilog = "Experimental version. Use at own risk." )
parser.add_option( "-d", "--distance", action="store", dest="radius",
default=0.015, metavar="DISTANCE", type="float",
help="displacement distanse [default: %default]" )
parser.add_option( "-r", "--rotate", action="store_true", dest="rotate",
default=False, help="place points horizontaly, if count = 2 [default: %default]" )
parser.add_option( "-l", "--log", action="store", dest="logFile",
metavar="FILE",
help="log overlapped points to file [default: %default]. Makes no sense without '-f' option" )
parser.add_option( "-f", "--field", action="store", dest="logField",
metavar="FIELD",
help="Shape field that will be logged to file. Makes no sense without '-l' option" )
( opts, args ) = parser.parse_args()
if len( args ) != 2:
parser.error( "incorrect number of arguments" )
if opts.logFile and opts.logField:
doDisplacement( args[ 0 ], args[ 1 ], opts.radius, opts.rotate, opts.logFile, opts.logField )
else:
doDisplacement( args[ 0 ], args[ 1 ], opts.radius, opts.rotate )
The script
- finds displacement candidates by drawing a circle around each point
- calculates the angle depending on the number of points to displace (e.g. 2 - 0 & 180)
- moves the points
I'm using this input file. Executing python displacement.py input.geojson output.geojson -d 100 should displace all points, but there are always 0 points displaced.
May you please help me finding out the cause?
