葵花8、9号卫星数据处理及可视化

总想着躺赢

已于 2024-06-19 13:56:44 修改

阅读量1.5k

点赞数 16

文章标签： python

于 2024-06-19 11:28:47 首次发布

本文链接：https://blog.csdn.net/qq_43352337/article/details/139792080

版权

H8/AHI数据读取及可视化

前一篇文章讲解了葵花8号卫星数据的预处理，考虑到属于遥感的内容，因此将数据处理为tiff格式。然而，在实际应用中发现，很多不会使用ENVI的小伙伴感觉tiff格式不好用。因此，对数据处理过程给出了新的处理方式。

本文主要内容：

1.H8/H9数据读取； 2. 生成可视化的真彩色图像。

注意：本文使用的数据来自H8官网的NC文件

废话不多说，直接上代码！

# 导入相关的库
import os
import numpy as np
from xarray import open_dataset, DataArray
import matplotlib.pyplot as plt
from cartopy import crs
from scipy.interpolate import interp1d
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# 定义文件
filename = r'D:\Satellite_Data\H8DATA\NC_H08_20170311_0500_R21_FLDK.02401_02401.nc'

# 使用xarray读取数据集
dataset = open_dataset(filename)

# 定义研究区域范围
lonmin = 80
lonmax = 140
latmin = 0
latmax = 60

# 读取经纬度信息
lon = dataset['longitude'].loc[lonmin:lonmax]
lat = dataset['latitude'].loc[latmax:latmin]

# 读取反射率数据(R:0.64um, G:0.51um, B:0.47um)
B = dataset['albedo_01'].loc[latmax:latmin, lonmin:lonmax]
G = dataset['albedo_02'].loc[latmax:latmin, lonmin:lonmax]
R = dataset['albedo_03'].loc[latmax:latmin, lonmin:lonmax]

# 通过这种方式读取出来的数据可以直接可视化绘图
B.plot()
plt.show()


# 将三个波段合成为RGB真彩色图像

## 1.对数据拉伸
def stretch(data):
    data = (data - 0) / (1 - 0) * 255
    x = [0, 30, 60, 120, 190, 255]
    y = [0, 110, 160, 210, 240, 255]
    interp = interp1d(x, y, bounds_error=False, fill_value=255)

    return interp(data).astype(int)


B = stretch(B.to_numpy())
G = stretch(G.to_numpy())
R = stretch(R.to_numpy())

RGB = np.dstack((R, G, B)) / 255
Color = RGB.reshape(-1, 3)

## 2.可视化绘图
figure = plt.figure(figsize=(8, 8))
proj = crs.PlateCarree()
ax = figure.add_subplot(111, projection=proj)
ax.coastlines()
ax.pcolormesh(lon, lat, RGB, color=Color, shading="nearest", transform=crs.PlateCarree())
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(lonmin, lonmax, 10), crs=proj)
ax.set_yticks(np.arange(latmax, latmin, -10), crs=proj)
plt.tight_layout()
plt.show()