python - Python OGR:在源代码中寻找代码位移的坏点
问题描述
我试图将一个 python ogr 脚本改编为 Python3,该脚本将点彼此太靠近(最初来自http://wiki.gis-lab.info/w/Displacement_of_points_with_same_coordinates,_ShiftPoints/OGR for python2)。我还将 GDAL 驱动程序从 Shapefile 更改为 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 )
剧本
- 通过在每个点周围画一个圆圈来查找候选位移
- 根据要置换的点数计算角度(例如 2 - 0 和 180)
- 移动点
我正在使用这个输入文件。执行python displacement.py input.geojson output.geojson -d 100
应该置换所有点,但总有 0 个点被置换。
请你帮我找出原因吗?
解决方案
这段代码没有坏处,它完全按照它所说的那样做。
它只置换具有完全相同坐标的点。您的文件有不同的坐标。我用具有相同坐标的 shapefile 进行了尝试,效果很好。
我想你将不得不寻找另一种解决方案。
推荐阅读
- macos - 如何修复在 MacOS 中从终端打开 sublime?
- javascript - 如何在功能组件中从父级更改子状态组件
- spring - 无法在 Spring Boot 应用程序中运行 Embedded Keycloak
- html - 反应解码字符串中的特殊字符 - 是否有任何内置机制或功能
- javascript - 变量/参数未定义
- c# - .Net Core 类库 - 复制文件以发布 (nupkg)
- angular - 如果 'router-outlet' 是一个 Angular 组件,那么验证它是这个模块的一部分
- flutter - Flutter TextFormField 禁用键盘但允许接受输入
- c++ - 为什么破坏顺序与构造顺序相同,使用静态对象(C++)?
- html - 如何模糊除悬停在 CSS 和 HTML 中的元素之外的所有内容