Python将GIS栅格转为CSV并网格绘点异常问题求助

拥有EPSG:7416坐标系的丹麦建筑栅格文件denmark_buildings.tif，通过rioxarray读取并绘制栅格图显示正常，但通过GDAL将其转为XYZ格式再保存为CSV后，读取坐标绘制散点图时，全区域布满红点，无法识别建筑形态。

操作步骤：

1. 读取栅格文件：

import rioxarray
import matplotlib.pyplot as plt
from osgeo import gdal, ogr, osr
import xarray
import rioxarray
import numpy as np
import rasterio

raster_reprojected = rioxarray.open_rasterio("denmark_buildings.tif").squeeze()

2. GDAL转XYZ格式：

# Convert using GDAL Translate Driver
# read .tif file
ds_buildings = gdal.Open("denmark_buildings.tif")
ds_buildings_xyz = gdal.Translate("denmark_buildings.xyz", ds_buildings)
ds_buildings_xyz = None

3. 读取CSV并绘制散点图：

## Read csv as numpy array
buildings_array = np.loadtxt("denmark_buildings.csv", delimiter=',', skiprows=1)
# Get X and Y coordinates as lists
x_buildings = buildings_array[:, 0]
y_buildings = buildings_array[:, 1]

# Create a figure and a set of subplots
fig, ax = plt.subplots(1, figsize=(8, 12), dpi=500)

# Plot the data
ax.scatter(x_buildings, y_buildings, s=0.01, c='red', label='Buildings')

min_longitude = 620561.9971020991
max_longitude = 621791.9681244957
min_latitude = 6108663.515810543
max_latitude = 6110367.123382721
# Zoom in to the area of interest
# Set the limits of x and y to the region you want to zoom into
ax.set_xlim([min_longitude, max_longitude])
ax.set_ylim([min_latitude, max_latitude])

# Show the plot
plt.show()

GDAL的Translate工具默认会输出栅格中所有像素的坐标与值，包括代表无建筑区域的背景（NoData）像素。散点图将所有这些像素点都绘制出来，自然会布满整个区域，掩盖了建筑对应的有效像素。

方案1：转换XYZ时过滤背景像素 #

先获取栅格的NoData值，转换时忽略这些背景像素：

ds_buildings = gdal.Open("denmark_buildings.tif")
# 获取栅格的NoData值
nodata_value = ds_buildings.GetRasterBand(1).GetNoDataValue()
# 转换时排除NoData像素
ds_buildings_xyz = gdal.Translate(
    "denmark_buildings.xyz", 
    ds_buildings, 
    srcNodata=nodata_value,
    dstNodata=nodata_value
)
ds_buildings_xyz = None

之后再将XYZ转为CSV，读取后绘制的散点图就只会包含建筑区域的像素点。

方案2：直接用rioxarray提取有效坐标（更高效） #

跳过GDAL转换步骤，直接从栅格中筛选出建筑对应的像素坐标：

import rioxarray
import numpy as np
import matplotlib.pyplot as plt

# 读取栅格并压缩维度
raster = rioxarray.open_rasterio("denmark_buildings.tif").squeeze()
# 生成掩码：筛选出非NoData的像素（即建筑区域）
building_mask = raster != raster.rio.nodata
# 获取对应像素的X、Y坐标
x_coords = raster.x[building_mask].values
y_coords = raster.y[building_mask].values

# 绘制散点图
fig, ax = plt.subplots(1, figsize=(8, 12), dpi=500)
ax.scatter(x_coords, y_coords, s=0.01, c='red', label='Buildings')
# 设置感兴趣区域范围
ax.set_xlim([620561.9971020991, 621791.9681244957])
ax.set_ylim([6108663.515810543, 6110367.123382721])
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者errenmike1806