python使用gdal将栅格文件转为shp
将栅格数据转为shp面
from osgeo import gdalconst, gdal, ogr, osr
import os
def raster2poly(raster, outshp):
inraster = gdal.Open(raster)
inband = inraster.GetRasterBand(1)
prj = osr.SpatialReference()
prj.ImportFromWkt(inraster.GetProjection())
drv = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outshp):
drv.DeleteDataSource(outshp)
Polygon = drv.CreateDataSource(outshp)
Poly_layer = Polygon.CreateLayer(raster[:-4], srs=prj, geom_type=ogr.wkbMultiPolygon)
newField = ogr.FieldDefn('value', ogr.OFTReal)
Poly_layer.CreateField(newField)
gdal.FPolygonize(inband, None, Poly_layer, 0)
Polygon.SyncToDisk()
Polygon = None
raster2poly("input_tif.tif", "output_shp.shp")
将栅格数据转为shp点
from osgeo import ogr, osr, gdal
import numpy as np
import pandas as pd
inputTIF = "input_tif.tif"
NoData_value = -999
result_path = 'output_shp.shp'
def createPoint_strNosrs(reader, result_path):
"""
创建无空间参考且字段类型全部为string的点shp文件
:param reader: 待生成shp文件的表属性信息 为dict格式 其中最后两个key为点的位置信息 key='point'
:param result_path: 待生成的shp 文件的存储地址
:return:
"""
driver = ogr.GetDriverByName("ESRI Shapefile")
data_source = driver.CreateDataSource(result_path)
layer = data_source.CreateLayer("point", geom_type=ogr.wkbPoint)
l_list = []
l_col = list(reader[0].keys())
for i in range(len(l_col) - 1):
l_list.append({list(reader[0].keys())[i]: ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString)})
ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString).SetWidth(120)
layer.CreateField(ogr.FieldDefn(str(l_col[i][:10]), ogr.OFTString))
for row in reader:
feature = ogr.Feature(layer.GetLayerDefn())
for use_i in l_list:
feature.SetField(str(list(use_i.keys())[0][:10]), row[list(use_i.keys())[0]])
wkt = row['point']
point = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(point)
layer.CreateFeature(feature)
feature = None
data_source = None
data = gdal.Open(inputTIF, gdal.GA_ReadOnly)
data1= data.GetRasterBand(1).ReadAsArray()
transform = data.GetGeoTransform()
res = np.where(data1 != NoData_value )
lon = transform[3] + res[1] * transform[4] + res[0] * transform[5]
lat = transform[0] + res[1] * transform[1] + res[0] * transform[2]
values = data1[data1 != NoData_value ]
point = ['POINT (' + str(lat[k]) + ' ' + str(lon[k])+')' for k in range(len(lon))]
dataUse = pd.DataFrame({"value": values, "point": point})
reader = dataUse.to_dict(orient='records')
createPoint_strNosrs(reader, result_path)