使用gdal_polygonize.py将TIFF文件转换为成型时，无法提取所有点

Question

我试图通过使用gdal_polygonize.py命令将其转换为linux上的

conversion成功，从而从TIFF文件中提取高程点，但是形状文件不包含所有高程点。

我正在使用以下命令进行转换

gdal_polygonize.py NT60ne_DTM_2m.tif -f "ESRI Shapefile" NT60ne_DTM_2m.shp -fieldname elevation

的屏幕截图。

以下是 nt60ne_dtm_2m.2m.tif 和 nt60ne_dtm_2m.sshp

nt60ne_dtm_2m.shp

我知道转换后的文件是不完整的，因为我使用 raster pixel在qGIS工具中执行了相同的操作及以下是其输出

使用 gdal_polygonize.py.py 命令进行转换时，我在这里缺少什么。为什么不完整？

update ：添加源文件供其他用户尝试 nt60ne_dtm_2m.tif

Answer 1

您要比较的两个功能不会产生相同的输出。 gdal_polygonize将创建连接等值像素（在单个多边形）的多边形。与QGIS函数相比，该功能将每个像素转换为唯一的点/多边形，无论其值如何。

我尚不清楚您的GDAL_Polygonize的结果似乎包含点而不是多边形。那应该是不可能的。

您可以尝试使用 gdal.polygonize 从Python进行同样的操作。当我使用下面的片段时，它似乎可以正常工作。输出是下图中的红线，在输入栅格上覆盖。

from osgeo import gdal, ogr, osr

ras_file = 'NT60ne_DTM_2m.tif'
vec_file = 'NT60ne_DTM_2m.gpkg'

ds_ras = gdal.OpenEx(ras_file)
src_band = ds_ras.GetRasterBand(1)
srs_wkt = ds_ras.GetProjection()

srs = osr.SpatialReference()
srs.ImportFromWkt(srs_wkt)


drv = ogr.GetDriverByName('GPKG')
ds_vec = drv.CreateDataSource(vec_file)
lyr = ds_vec.CreateLayer('', srs=srs, geom_type=ogr.wkbPolygon)

lyr.CreateField(ogr.FieldDefn("elevation", ogr.OFTReal))
fld = lyr.GetLayerDefn().GetFieldIndex("elevation")

gdal.Polygonize(src_band, None, lyr, fld)


srs = None
ds_ras = None
ds_vec = None

一个子集，仅显示下右角：

编辑：

下面的摘要是将每个像素写入功能（point/polygon）到磁盘上的示例，并将栅格值附加为属性。您可以避免使用Shapely＆amp;如果需要的话，地理器需要更多的代码才能仅使用OGR来创建向量文件。

可以使用相同的技术来输出多边形，但创建

这样的方法，将字符串格式化/解析用于大量点永远不会表现出色。

import geopandas as gpd
from shapely import wkt
import pandas as pd
import numpy as np
from osgeo import gdal, ogr, osr

ras_file = 'NT60ne_DTM_2m.tif'
vec_file = 'NT60ne_DTM_2m_points.gpkg'

# read raster data & metadata
ds = gdal.OpenEx(ras_file)
gt = ds.GetGeoTransform()
xs = ds.RasterXSize
ys = ds.RasterYSize
srs_wkt = ds.GetProjection()
nodata = ds.GetRasterBand(1).GetNoDataValue()
data = ds.ReadAsArray()
ds = None

# mask to ignore nodata
assert np.isfinite(nodata) # else use np.isfinite / np.isnan
valid_data = data != nodata

# (outer) extent of raster
ulx, xres, _, uly, _, yres = gt
lrx = ulx + xs * xres
lry = uly + ys * yres

# convert extent to center of pixel
ulx = ulx + xres/2
lrx = lrx - xres/2
uly = uly + yres/2
lry = lry - yres/2

# create coordinate grids
mapy, mapx = np.mgrid[uly:lry+yres:yres, ulx:lrx+xres:xres]

# create normal DataFrame
df = pd.DataFrame(dict(mapx=mapx[valid_data], mapy=mapy[valid_data], elevation=data[valid_data]))

# Use WKT Point to add the geometry (just a string column)
df["geometry"] = df.apply(lambda row: wkt.loads(f"POINT ({row.mapx} {row.mapy})"), axis=1)

# convert to actual GeoDataFrame
df = gpd.GeoDataFrame(df, geometry='geometry', crs=srs_wkt)

# save to disk
df.to_file(vec_file, driver="GPKG")