使用 matplotlib 底图绘制 GDAL 栅格
Posted
技术标签:
【中文标题】使用 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。以上是关于使用 matplotlib 底图绘制 GDAL 栅格的主要内容,如果未能解决你的问题,请参考以下文章
没有河流的世界地图与matplotlib / Basemap?
GeoViews:使用 matplotlib 后端添加平铺底图