python浮雕图片_Python+basemap绘制投影图像,PythonBasemap,已,的,影像

最新推荐文章于 2023-07-19 08:00:00 发布

weixin_39946429

最新推荐文章于 2023-07-19 08:00:00 发布

阅读量237

点赞数

文章标签： python浮雕图片

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)