+关注
已关注

分类  

暂无分类

标签  

暂无标签

日期归档  

2019-07(1)

2019-08(109)

2019-09(120)

2019-10(17)

2019-11(1)

Python+Basemap绘制已投影的影像

发布于2020-09-28 19:46     阅读(909)     评论(0)     点赞(26)     收藏(2)


0

1

2

3

4

5

前言

(1) 数据是一副具有albers投影的地温影像,而Basemap的基本参数中要求输入影像的左下右上的经纬度坐标,所以关键在于如何将投影坐标转化为大地经纬度坐标。
(2) 在坐标转换过程中,用到了pyproj库,首先要定义转换前后的坐标类型,定义WGS84坐标系可用pyproj.CRS.from_epsg(4326)进行定义,但是由于我的pyproj的数据库路径有些问题,所以首先采用gdal读取一幅WGS84坐标的影像的wkt投影信息,然后pyproj.CRS.from_wkt( )进行定义要转换成的坐标系。
(3) 关于影像原有的albers投影信息,我是采用Proj4定义的,如果你的原始影像的投影信息不知如何用proj4进行表达,可以用我定义WGS84坐标信息的方法来定义,或者用gdal读取投影信息后输出proj4信息。更多关于gdal读取输出影像投影信息的内容可以参考我之前的博客Python GDAL学习笔记(二)
(4) 为了更有实用性,绘制了多(两)个子图。
(5) 代码中涉及仿射变换参数的知识,可参考以下。

dfGeoTransform[0] /* top left x 左上角x坐标*/
adfGeoTransform[1] /* w--e pixel resolution 东西方向上的像素分辨率*/
adfGeoTransform[2] /* rotation, 0 if image is "north up" 如果北边朝上,地图的旋转角度*/
adfGeoTransform[3] /* top left y 左上角y坐标*/
adfGeoTransform[4] /* rotation, 0 if image is "north up" 如果北边朝上,地图的旋转角度*/
adfGeoTransform[5] /* n-s pixel resolution 南北方向上的像素分辨率*/

代码

from osgeo import gdal
from osgeo import gdal_array
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# import rasterio
import pyproj
import math

file = r'F:\***.tif'

# 读取tif影像函数            
def readTifAsArray(tifPath):
    dataset = gdal.Open(tifPath)     
    if dataset == None:
        print(tifPath + "文件错误")
        return tifPath
            
    image_datatype = dataset.GetRasterBand(1).DataType  
    row = dataset.RasterYSize
    col = dataset.RasterXSize
    nb  = dataset.RasterCount
    proj = dataset.GetProjection()
    gt = dataset.GetGeoTransform()
    
    if nb != 1:
        array = np.zeros((row, col, nb), 
                        dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
                                image_datatype))
        for b in range(nb):
            band = dataset.GetRasterBand(b + 1)
            nan = band.GetNoDataValue()
            array[:, :, b] = band.ReadAsArray()
    else:
        array = np.zeros((row,col),
                         dtype = gdal_array.GDALTypeCodeToNumericTypeCode(
                         image_datatype))
        band = dataset.GetRasterBand(1)        
        nan = band.GetNoDataValue()        
        array = band.ReadAsArray()
    return array, nan, gt, proj

data  = readTifAsArray(file)[0] # 栅格数组
data[data>320] = np.nan # 去除海面 999
affine = readTifAsArray(file)[2] # 读取仿射变换参数
width = data.shape[1] 
height = data.shape[0]

map_width = width * affine[1] # 影像的宽度米
map_height = height * affine[1] # 高度 米
xmin = affine[0] # 分别为左下xy 右上xy坐标
xmax = xmin + map_width
ymax = affine[3]
ymin = ymax - map_height
extent = [xmin, xmax, ymin, ymax] # [left, right, bottom, top]

crs = readTifAsArray(file)[3] # 原始影像的投影albers
crs2 = readTifAsArray(r'Y:\folder\仿射变换_final_test.tif')[3] # WGS84大地坐标系

crs = pyproj.CRS.from_proj4('+proj=aea +lat_0=0 +lon_0=105 +lat_1=25 +lat_2=47 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs')
crs2 = pyproj.CRS.from_wkt(crs2)

x1 = np.array([xmin, xmax]) # 一定要转换为numpy数组形式,否则转换后可能出现inf等无效值
y1 = np.array([ymin, ymax])

proj = pyproj.Transformer.from_crs(crs, crs2)
lat, lon = proj.transform(x1, y1) # 将左下右上xy坐标(米)转化为大地坐标系(经纬度)

# Instantiate Basemap class
plt.style.use(['no-latex'])
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
fig = plt.figure(figsize=(8,5))
ax = fig.add_subplot(121)

m = Basemap(llcrnrlon=lon[0], llcrnrlat=lat[0], urcrnrlon=lon[1], urcrnrlat=lat[1], # 左下右上经纬度坐标 关键在于如何将已投影的坐标转换为经纬度坐标
            projection='aea', # albers等面积投影
            resolution='h', lat_0=0, lon_0=105, lat_1=25, lat_2=47) # 指定albers的参考纬度线和中心经度线
            # There might be other parameters to set depending on your CRS

# m.drawmapboundary(fill_color='skyblue') # 背景绘制蓝色
# m.shadedrelief() # 灰蓝背景
# m.etopo() # 绘制noaa浮雕底图 etopo
# m.bluemarble() # 大理石蓝背景
# m.fillcontinents(lake_color='aqua')
# m.drawcoastlines() # 绘制世界矢量边界

import matplotlib
norm = matplotlib.colors.Normalize(vmin=275, vmax=325)
m.imshow(data, origin='upper', extent=extent, cmap='jet', norm=norm) # 绘制栅格数据

cb = m.colorbar(extend='both' ,location="bottom", pad=0.3) 

cb.ax.tick_params(labelsize=10)  #设置色标刻度字体大小。
font = {'family' : 'serif',
#       'color'  : 'darkred',
    'color'  : 'black',
    'weight' : 'normal',
    'size'   : 10,
    }
cb.set_label('Temperature/K' ,fontdict=font)

parallels = np.arange(0.,43,1.) # 绘制经纬度线
meridians = np.arange(116.,125.,2.)

m.drawparallels(parallels,labels=[True,False,False,False], linewidth=0.1) # labels = [left,right,top,bottom]
m.drawmeridians(meridians,labels=[False,False,False,True], linewidth=0.1)
# plt.show()
# plt.savefig('figure_name.png',dpi=300, bbox_inches='tight') # , transparent=True透明保存

############################################################################### 绘制第二个子图

ax = fig.add_subplot(122)
m = Basemap(llcrnrlon=lon[0], llcrnrlat=lat[0], urcrnrlon=lon[1], urcrnrlat=lat[1], # 左下右上经纬度坐标 关键在于如何将已投影的坐标转换为经纬度坐标
            projection='aea', # albers等面积投影
            resolution='h', lat_0=0, lon_0=105, lat_1=25, lat_2=47) # 指定albers的参考纬度线和中心经度线

import matplotlib
norm = matplotlib.colors.Normalize(vmin=275, vmax=325)
m.imshow(data, origin='upper', extent=extent, cmap='jet', norm=norm) # 绘制栅格数据

cb = m.colorbar(extend='both' ,location="bottom", pad=0.3) 

cb.ax.tick_params(labelsize=10)  #设置色标刻度字体大小。
font = {'family' : 'serif',
#       'color'  : 'darkred',
    'color'  : 'black',
    'weight' : 'normal',
    'size'   : 10,
    }
cb.set_label('Temperature/K' ,fontdict=font)

parallels = np.arange(0.,43,1.) # 绘制经纬度线
meridians = np.arange(116.,125.,2.)

m.drawparallels(parallels,labels=[True,False,False,False], linewidth=0.1) # labels = [left,right,top,bottom]
m.drawmeridians(meridians,labels=[False,False,False,True], linewidth=0.1)
plt.savefig('figure_name.png', dpi=300, bbox_inches='tight') # , transparent=True透明保存

结果

在这里插入图片描述

原文链接:https://blog.csdn.net/tk20190411/article/details/108821478

0

1

2

3

4



所属网站分类: 技术文章 > 博客

作者:你太美丽

链接: https://www.pythonheidong.com/blog/article/550864/1143e6f96a1445796013/

来源: python黑洞网

任何形式的转载都请注明出处,如有侵权 一经发现 必将追究其法律责任

26 0
收藏该文
已收藏

评论内容:(最多支持255个字符)