目录
一.有些遥感影像虽然是tif格式的,但是直接运行会报错,把遥感影像在ArcGIS中重新导出数据,适用于国产数据
1、不要在anaconda中直接安装GDAL,先下一个XXX.whl在安装,参考上面的链接
3、用python打包成exe的时候,python 3.9 + gdal 3.3.2版本会报错,回退到python3.7 + gdal 3.3.1就好了(提醒打包失败,一般都是环境的问题)
随机森林代码参考,在此表示感谢
代码运行遇到的问题:
一.有些遥感影像虽然是tif格式的,但是直接运行会报错,把遥感影像在ArcGIS中重新导出数据,适用于国产数据
二.osgeo的gdal包的安装问题:
1、不要在anaconda中直接安装GDAL,先下一个XXX.whl在安装,参考上面的链接
2、Cannot find proj.db的错误
ERROR 1: PROJ: proj_identify: Cannot find proj.db
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db
解决方法一:把GDAL中的proj文件添加到系统变量中
解决方法二:
如果打包成exe换个电脑,就又会缺少环境变量,又会报这个错误,
解决办法,把proj文件复制到exe同一文件下,在代码中添加如下代码:
path1 = sys.executable
path2 = os.path.dirname(path1)
path = os.path.join(path2, 'proj')
os.environ['PROJ_LIB'] = path
print(os.environ['PROJ_LIB'])
3、用python打包成exe的时候,python 3.9 + gdal 3.3.2版本会报错,回退到python3.7 + gdal 3.3.1就好了(提醒打包失败,一般都是环境的问题)
Rater转MultiShp
小图斑去除
其实做一个闭运算更好,但是GDAL没提供具体的算法,只能用cv2的算法,请自行百度把
from osgeo import gdal
def callback(v1, v2, v3):
print(v1)
# Remove small spots
def RemoveSpots(inRaster, threshold=20, connectedness=4):
print(inRaster)
inputFile = inRaster
options = ()
mask = 'none'
src_ds = gdal.Open(inputFile, gdal.GA_Update)
srcBand = src_ds.GetRasterBand(1)
dstBand = srcBand
if mask is 'default':
maskBand = srcBand.GetMaskBand()
elif mask is 'none':
maskBand = None
else:
mask_ds = gdal.Open(mask)
maskBand = mask_ds.GetRasterBand(1)
# 参数说明 输入数据波段 、设置掩码波段(只对掩码区域进行处理)、输出数据波段、去除板块的最大像元个数、图斑连通方式、占位(options)、进度条
result = gdal.SieveFilter(srcBand, maskBand, dstBand, threshold, connectedness, options, callback)
return result
if __name__ == '__main__':
gdal.AllRegister()
inRaster = "C:/Users/z6q6k6/Desktop/aa/result_1.tif"
RemoveSpots(inRaster)
栅格矢量化
from osgeo import gdal
from osgeo import ogr
from pathlib import Path
from osgeo import osr
import os
import sys
def callback(v1, v2, v3):
print(v1)
def RasterToShp(inputFile, dst_filename):
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(dst_filename):
driver.DeleteDataSource(dst_filename)
ds = gdal.Open(inputFile, gdal.GA_ReadOnly)
srcband = ds.GetRasterBand(1)
maskband = srcband.GetMaskBand()
prj = osr.SpatialReference()
prj.ImportFromWkt(ds.GetProjection())
drv = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = drv.CreateDataSource(dst_filename)
dst_layername = 'out'
dst_layer = dst_ds.CreateLayer(dst_layername, prj, ogr.wkbMultiPolygon)
dst_fieldname = 'Value'
fd = ogr.FieldDefn(dst_fieldname, ogr.OFTInteger)
dst_layer.CreateField(fd)
dst_field = 0
options = []
# 参数 输入栅格图像波段\掩码图像波段、矢量化后的矢量图层、需要将DN值写入矢量字段的索引、算法选项、进度条回调函数、进度条参数
print("栅格转shp完成")
return gdal.Polygonize(srcband, maskband, dst_layer, dst_field, options, callback)
if __name__ == '__main__':
gdal.AllRegister()
inputFile = "C:/Users/z6q6k6/Desktop/aa/result_1.tif"
dst_filename = "C:/Users/z6q6k6/Desktop/aa/allShp.shp"
RasterToShp(inputFile, dst_filename)
按分类结果提取每一类的矢量文件
把包含的所有类别的shp复制多份,然后在每一份中,删除不符合要求的类别,就得到了每一类的矢量文件
这个方法目前看起来比较low,但是还是学到了很多东西,值得分享!
在这再提供一种思路:
getValue函数中,feature.GetFieldAsString('Value')获取类别
然后参考CopyShp函数,分别导入到新创建的矢量中
from osgeo import gdal
from osgeo import ogr
from pathlib import Path
from osgeo import osr
import os
import sys
import shutil
def CopyShp(inputfile, pathfile, str):
os.chdir(pathfile)
# 设置driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# 打开输入的矢量
inDs = driver.Open(inputfile, 1)
if inDs is None:
print("Could not open", inputfile)
sys.exit(1)
# 打开输入矢量的图册
inLayer = inDs.GetLayer()
# 检查所需创建的矢量是否已存在
if os.path.exists(str):
driver.DeleteDataSource(str)
# 创建矢量
outDs = driver.CreateDataSource(str)
if outDs is None:
print("Could not create file")
sys.exit(1)
# 创建图册
prj = osr.SpatialReference()
a = inLayer.GetSpatialRef()
s = a.ExportToWkt()
prj.ImportFromWkt(s)
outLayer = outDs.CreateLayer(str, srs=prj, geom_type=ogr.wkbPolygon)
# for inFeature in inLayer:
# outLayer.CreateFeature(inFeature)
# 将输入矢量的属性应用到新矢量
fieldDefn = inLayer.GetFeature(0).GetFieldDefnRef('Value')
outLayer.CreateField(fieldDefn)
featureDefn = outLayer.GetLayerDefn()
inFeature = inLayer.GetNextFeature()
while inFeature:
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(inFeature.GetGeometryRef())
for i in range(inFeature.GetFieldCount()):
value = inFeature.GetField(i)
outFeature.SetField(i, value)
outLayer.CreateFeature(outFeature)
inFeature.Destroy()
outFeature.Destroy()
inFeature = inLayer.GetNextFeature()
inDs.Destroy()
outDs.Destroy()
return os.path.join(pathfile, str)
def repack(shpFileName, strValue):
driver = ogr.GetDriverByName('ESRI Shapefile')
pFeatureDataset = driver.Open(shpFileName, 1)
if pFeatureDataset is None:
print("Could not open", pFeatureDataset)
sys.exit(1)
pFeaturelayer = pFeatureDataset.GetLayer()
# strValue = 0
strFilter = "Value != '" + str(strValue) + "'"
pFeaturelayer.SetAttributeFilter(strFilter)
pFeatureDef = pFeaturelayer.GetLayerDefn()
pLayerName = pFeaturelayer.GetName()
pFieldName = "Value"
# pFieldIndex = pFeatureDef.GetFieldIndex(pFieldName)
pFeatureDef.GetFieldIndex(pFieldName)
for pFeature in pFeaturelayer:
pFeatureFID = pFeature.GetFID()
pFeaturelayer.DeleteFeature(int(pFeatureFID))
strSQL = "REPACK " + str(pFeaturelayer.GetName())
pFeatureDataset.ExecuteSQL(strSQL, None, "")
pFeatureLayer = None
pFeatureDataset = None
return 0
def getValue(inShp):
gdal.SetConfigOption("SHAPE_ENCODING", "")
ogr.UseExceptions()
in_ds = ogr.Open(inShp)
in_lyr = in_ds.GetLayer()
list = []
for feature in in_lyr:
# print(feature.SetAttributeFilter("Value = '1'"))
if not feature.GetFieldAsString('Value') in list:
list.append(feature.GetFieldAsString('Value'))
del in_ds
return list
def SimpleToMulti(inputfile, pathfile):
if os.path.exists(pathfile):
# os.remove(outfile) 删除的时非空文件夹,os.remove会出现拒绝访问
shutil.rmtree(pathfile)
os.makedirs(pathfile)
list = getValue(inputfile)
sum = len(list)
num = 0
for i in list:
pathAndName = CopyShp(inputfile, pathfile, i + '.shp')
repack(pathAndName, i)
num = num + 1
print(f"处理完成{round(num / sum, 2) * 100}%")
if __name__ == '__main__':
gdal.AllRegister()
inputfile = "C:/Users/z6q6k6/Desktop/aa/allShp.shp"
pathfile = "C:/Users/z6q6k6/Desktop/aa//eachShp"
SimpleToMulti(inputfile, pathfile)
矢量平滑
用两个缓冲区(正缓冲区和负缓冲区)来平滑矢量的线
但是有一个问题,就是如果geometry.GetArea()太大,速度会很慢
from osgeo import gdal
from osgeo import ogr
from pathlib import Path
import tqdm
from osgeo import osr
import os
import sys
def callback(v1, v2, v3):
print(v1)
def smoothing(inShp, fname, bdistance=0.001):
"""
:param inShp: 输入的矢量路径
:param fname: 输出的矢量路径
:param bdistance: 缓冲区距离
:return:
"""
ogr.UseExceptions()
in_ds = ogr.Open(inShp)
in_lyr = in_ds.GetLayer()
# 创建输出Buffer文件
driver = ogr.GetDriverByName('ESRI Shapefile')
if Path(fname).exists():
driver.DeleteDataSource(fname)
# 新建DataSource,Layer
out_ds = driver.CreateDataSource(fname)
out_lyr = out_ds.CreateLayer(fname, in_lyr.GetSpatialRef(), ogr.wkbPolygon)
def_feature = out_lyr.GetLayerDefn()
# 遍历原始的Shapefile文件给每个Geometry做Buffer操作
for feature in tqdm.tqdm(in_lyr):
geometry = feature.GetGeometryRef()
# print(geometry.GetArea())
buffer = geometry.Buffer(bdistance).Buffer(-bdistance)
out_feature = ogr.Feature(def_feature)
out_feature.SetGeometry(buffer)
out_lyr.CreateFeature(out_feature)
out_feature = None
out_ds.FlushCache()
del in_ds, out_ds
if __name__ == '__main__':
gdal.AllRegister()
inputshp = "C:/Users/z6q6k6/Desktop/aa/eachShp/2.shp"
output_path = "C:/Users/z6q6k6/Desktop/aa/eachShp/2_1.shp"
smoothing(inputshp, output_path)
print("end!")
异步问题!!!
在读取文件时候多是异步操作,如果上一步的文件还未完全创建(可能只有几k),就立即进行下一步的操作,一定会报错
感谢
本文参考了很多博客,不能一一添加引用,在此一并表示感谢