如何使用 python 在特定区域 (lat/lon) 绘制 geotiff 数据

Question

我有一个带有高程数据初始化的geotiff栅格数据集，我想在特定区域绘制它，例如60°E - 70° E，70°S - 80°E。

我有一些来自 here,but the pcolormesh seem couldn't plot my geotif.it's all red. picture. The picture is shown by imshow as really picture

的代码

当我尝试使用下面的代码绘图时：

path = "F:\Mosaic_h1112v28_ps.tif"
dataset = gdal.Open(path)
data = dataset.ReadAsArray()
x0, dx, dxdy, y0, dydx, dy = dataset.GetGeoTransform()
nrows, ncols = data.shape
londata = np.linspace(x0, x0+dx*ncols)
latdata = np.linspace(y0, y0+dy*nrows)
lons, lats = np.meshgrid(lonarray, latarray)  
fig = plt.figure(figsize=(8, 8))
m = Basemap(projection='lcc', lon_0=67.5, lat_0=-68.5, height=950000,
            width=580000, resolution='h') 
m.drawcoastlines() 
x, y = m(lons, lats)

那我就不知道怎么继续了。我只想使用 imshow，但是 imshow 没有指定区域（lat/lon）。

非常感谢您的帮助。

Answer 1

这是我的解决方案。

1。导入 GEOTIFF 文件并将其转换为二维数组数据

from osgeo import gdal
pathToRaster = r'./xxxx.tif'
raster = gdal.Open(pathToRaster,  gdal.GA_ReadOnly)
data = raster.GetRasterBand(1).ReadAsArray()
data = data[::-1]

2。使用 Pcolormesh

绘制

kk = plt.pcolormesh(data,cmap = plt.cm.Reds,alpha = 0.45, zorder =2)

Answer 2

这是一个很好的问题，这是我的解决方案。

所需的软件包：georaster 及其依赖项（gdal 等）。
可从 http://dwtkns.com/srtm/

下载用于演示目的的数据

import georaster
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(8,8))

# full path to the geotiff file
fpath = r"C:\path_to_your\geotiff_file\srtm_57_10.tif"  # Thailand east

# read extent of image without loading
# good for values in degrees lat/long
# geotiff may use other coordinates and projection
my_image = georaster.SingleBandRaster(fpath, load_data=False)

# grab limits of image's extent
minx, maxx, miny, maxy = my_image.extent

# set Basemap with slightly larger extents
# set resolution at intermediate level "i"
m = Basemap( projection='cyl', \
            llcrnrlon=minx-2, \
            llcrnrlat=miny-2, \
            urcrnrlon=maxx+2, \
            urcrnrlat=maxy+2, \
            resolution='i')

m.drawcoastlines(color="gray")
m.fillcontinents(color='beige')

# load the geotiff image, assign it a variable
image = georaster.SingleBandRaster( fpath, \
                        load_data=(minx, maxx, miny, maxy), \
                        latlon=True)

# plot the image on matplotlib active axes
# set zorder to put the image on top of coastlines and continent areas
# set alpha to let the hidden graphics show through
plt.imshow(image.r, extent=(minx, maxx, miny, maxy), zorder=10, alpha=0.6)

plt.show()

结果图：

编辑1

我原来的回答着重于如何使用 Basemap 在最基本的投影上绘制简单的 geotiff 图像。如果无法访问所有必需的资源（即 geotiff 文件），就不可能有更好的答案。

这里我尝试改进我的答案。

我从整个世界的 geotiff 文件中剪下了一小部分。然后将其重新投影（变形）为 Basemap() 定义的 LCC 投影规范以供使用。所有的过程都是用 GDAL 软件完成的。生成的文件名为 "lcc_2.tiff"。使用这个 geotiff 文件，图像的绘制是用下面的代码完成的。

最重要的部分是geotiff文件必须与Basemap使用的投影具有相同的坐标系（相同的投影）。

import georaster
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(8,8))

m = Basemap(projection='lcc', lon_0=67.5, lat_0=-68.5, \
            height=950000, width=580000, resolution='h')

m.drawcoastlines()
m.fillcontinents(color='beige')

image = georaster.SingleBandRaster( "lcc_2.tiff", latlon=False)
plt.imshow(image.r, extent=image.extent, zorder=10, alpha=0.6)

plt.show()

输出图：

Answer 3

您可以使用 rioxarray

import rioxarray as rio
ds = rio.open_rasterio(path)
# Example lat lon range for subset
geometries = [
{
    'type': 'Polygon',
    'coordinates': [[
        [33.97301017665958, -118.45830810580743],
        [33.96660083660732, -118.37455741054782],
        [33.92304171545437, -118.37151348516299],
        [33.915042933806724, -118.42909440702563]
    ]]
}
]
clipped = ds.rio.clip(geometries)
clipped.plot()

如何使用 python 在特定区域 (lat/lon) 绘制 geotiff 数据

how to plot geotiff data in specific area (lat/lon) with python

python

gdal

matplotlib-basemap

1。导入 GEOTIFF 文件并将其转换为二维数组数据

2。使用 Pcolormesh