Gdal之栅格数据等高线

Posted 2023-04-13

参考技术A 1). 语法：

gdal_contour [ - b < band > ] [ - a < attribute_name > ] [ - amin < attribute_name > ] [ - amax < attribute_name > ]

             [ - 3d] [ - inodata]

             [ - snodata n] [ - i < interval > ]

             [ - f < formatname > ] [[ - dsco NAME = VALUE] ... ] [[ - lco NAME = VALUE] ... ]

             [ - off < offset > ] [ - fl < level >   < level >... ] [ - e < exp_base > ]

             [ - nln < outlayername > ] [ - q] [ - p]

              < src_filename >   < dst_filename >

2).示例 ：

案例：栅格数据生成等高线

场景： 想要给 d ata.tif 格式文件 生成10米等高线的 data.shp 格式文件

脚本：gdal_contour.exe -a elev data.tif data.shp -i 10.0

参考：https://www.osgeo.cn/gdal/programs/gdal_contour.html

脚本：BAT脚本示范

@echo off

chcp 936

@echo "gdal\gdal_contour.exe" -a elev data.tif data.shp -i 10.0

"gdal\gdal_contour.exe" -a elev data.tif data.shp -i 10.0

pause

使用 matplotlib 底图绘制 GDAL 栅格

【中文标题】使用 matplotlib 底图绘制 GDAL 栅格

【英文标题】：Plot GDAL raster using matplotlib Basemap

【发布时间】：2013-12-27 15:11:49

【问题描述】：

我想使用 matplotlib 底图绘制光栅 tiff (download-723Kb)。我的栅格投影坐标以米为单位：

In  [2]:
path = r'albers_5km.tif'
raster = gdal.Open(path, gdal.GA_ReadOnly)
array = raster.GetRasterBand(20).ReadAsArray()

print ('Raster Projection:\n', raster.GetProjection())
print ('Raster GeoTransform:\n', raster.GetGeoTransform())

Out [2]:
Raster Projection:
 PROJCS["unnamed",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",15],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",30],PARAMETER["longitude_of_center",95],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
Raster GeoTransform:
 (190425.8243, 5000.0, 0.0, 1500257.0112, 0.0, -5000.0)

如果我尝试使用 contourf 和 latlon=False 的 Robin 投影来绘制此图，则假定 x 和 y 是地图投影坐标（请参阅 docs，我认为这就是我所拥有的）。

但是如果我看一下情节，我会发现它的左下角非常小：

使用此代码：

In  [3]:
xy = raster.GetGeoTransform() 
x = raster.RasterXSize 
y = raster.RasterYSize    
lon_start = xy[0] 
lon_stop = x*xy[1]+xy[0] 
lon_step = xy[1]    
lat_start = xy[3] 
lat_stop = y*xy[5]+xy[3] 
lat_step = xy[5]

fig = plt.figure(figsize=(16,10)) 
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0)

lons = np.arange(lon_start, lon_stop, lon_step) 
lats = np.arange(lat_start, lat_stop, lat_step)    
xx, yy = np.meshgrid(lons,lats)

levels = [array.min(),-0.128305,array.max()] 
map.contourf(xx, yy,array, levels, cmap=cm.RdBu_r, latlon=False)

map.colorbar(cntr,location='right',pad='10%')    
map.drawcoastlines(linewidth=.5) 
map.drawcountries(color='red')

最终我不想有一个世界观，而是一个详细的观点。但这给了我一个绘制海岸线和国家的缩放级别，但数据再次放置在左下角，但不像上次那么小：

使用以下代码：

In  [4]:
extent = [ xy[0],xy[0]+x*xy[1], xy[3],xy[3]+y*xy[5]]
width_x = (extent[1]-extent[0])*10
height_y = (extent[2]-extent[3])*10

fig = plt.figure(figsize=(16,10))
map = Basemap(projection='stere', resolution='c', width = width_x , height = height_y, lat_0=40.2,lon_0=99.6,)

xx, yy = np.meshgrid(lons,lats)
levels = [array.min(),-0.128305,array.max()]
map.contourf(xx, yy, array, levels, cmap=cm.RdBu_r, latlon=False)

map.drawcoastlines(linewidth=.5)
map.drawcountries(color='red')

【问题讨论】：

不应该是 projection='stereo' 吗？ （缺少一个“o”可能只是一个错字）

不，它说“立体声”是不受支持的投影。 'stere' 用于立体投影

奇怪，如果我要缩写为 stereographic，它会是 stereo 而不是 stere ;-) 但也许所有投影的缩写都是 5 个字符左右。

您正在绘制 Albers (AEA) 坐标，就好像它们是 Robinson 或 Stereographic。在绘图之前，您需要将 xx 和 yy 重新投影到您的地图投影。

@RutgerKassies。感谢您指出这一点！我以为它是在即时执行此操作。我可以使用底图重新投影我的 xx 和 yy 吗？其次，当我将投影从“立体”更改为“aea”时，我希望它会正确显示，但仍如上所示显示。

【参考方案1】：

您可以使用以下代码转换坐标，它会自动将栅格的投影作为源，将 Basemap 对象的投影作为目标坐标系。

进口

from mpl_toolkits.basemap import Basemap
import osr, gdal
import matplotlib.pyplot as plt
import numpy as np

坐标转换

def convertXY(xy_source, inproj, outproj):
    # function to convert coordinates
    
    shape = xy_source[0,:,:].shape
    size = xy_source[0,:,:].size

    # the ct object takes and returns pairs of x,y, not 2d grids
    # so the the grid needs to be reshaped (flattened) and back.
    ct = osr.CoordinateTransformation(inproj, outproj)
    xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T))

    xx = xy_target[:,0].reshape(shape)
    yy = xy_target[:,1].reshape(shape)
    
    return xx, yy

读取和处理数据

# Read the data and metadata
ds = gdal.Open(r'albers_5km.tif')

data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

xres = gt[1]
yres = gt[5]

# get the edge coordinates and add half the resolution 
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5

ds = None

# create a grid of xy coordinates in the original projection
xy_source = np.mgrid[ymax+yres:ymin:yres, xmin:xmax+xres:xres]

绘图

# Create the figure and basemap object
fig = plt.figure(figsize=(12, 6))
m = Basemap(projection='robin', lon_0=0, resolution='c')

# Create the projection objects for the convertion
# original (Albers)
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)

# Get the target projection from the basemap object
outproj = osr.SpatialReference()
outproj.ImportFromProj4(m.proj4string)

# Convert from source projection to basemap projection
xx, yy = convertXY(xy_source, inproj, outproj)

# plot the data (first layer)
im1 = m.pcolormesh(xx, yy, data[0,:,:], cmap=plt.cm.jet, shading='auto')

# annotate
m.drawcountries()
m.drawcoastlines(linewidth=.5)

plt.savefig('world.png',dpi=75)

如果您需要像素位置 100% 正确，您可能需要自己检查坐标数组的创建非常小心（因为我根本没有）。希望这个例子能让你走上正轨。

【讨论】：

这确实让我走上了正轨。对于这个出色的答案，所有功劳都归功于您。能够在 Python 环境中进行分析和绘图真是太好了。

亲爱的 Rutger Kassies，这种方法在使用 Cartopy 时也有效吗？

如果光栅图像中只有 1 个波段，我可以知道上面代码中需要更改的内容吗？

@SaiKiran，这可能取决于你有多少维度。但如果你是 2D 的，你可能需要将这个：data[0,:,:] 改为 data。

@GabrielLucas，我​​对帖子进行了一些编辑。用np.mgrid翻转坐标生成，去掉数据的转置，加上shading="auto"。这消除了对我的警告。如果底图不适合您，请查看 Cartopy。