绘制出来的卫星云图

全圆盘真彩图:

全圆盘单通道红外图:

数据准备

数据格式说明:http://fy4.nsmc.org.cn/data/cn/data/realtime.html
数据下载地址:http://satellite.nsmc.org.cn/portalsite/Data/DataView.aspx?SatelliteType=1&SatelliteCode=FY4A#

本人使用的是4000M的全圆盘数据,下载数据需要申请账号

解析HDF数据

from netCDF4 import Dataset
import h5py# 两种解析方式 netCDF4 和 h5py   直接用conda安装hdf_data_path = "/Users/Downloads/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200805080000_20200805081459_4000M_V0001.HDF"# h5py 解析
hdf_obj = h5py.File(hdf_data_path, "r")
#打印文件里所有属性  属性含义自行查看数据说明格式
print(hdf_obj.keys())
# 通道1数据
print(hdf_obj['NOMChannel01'][:])# netCDF4 解析
nc_obj = Dataset(hdf_data_path)
#打印文件里所有属性
print(nc_obj.variables.keys())
# 通道1数据
print(nc_obj.variables['NOMChannel01'][:])

绘制单通道云图

from netCDF4 import Dataset
import matplotlib.pyplot as plthdf_data_path = "/Users/Downloads/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200805080000_20200805081459_4000M_V0001.HDF"nc_obj = Dataset(hdf_data_path)
type = nc_obj.variables.keys()
print(type)keyword = "NOMChannel"for k in type:if str(k).find(keyword) == 0:data = nc_obj.variables[k][:]plt.imshow(data, cmap='gray')plt.title(k)plt.axis('off')plt.show()print("通道" + k + "绘制完成")

运行绘制出14个通道图:

绘制真彩云图

官方卫星真彩云图:http://fy4.nsmc.org.cn/portal/cn/theme/FY4A.html

上面是官方绘制出来的,绘制这样一张云图需要多通道数据融合,搜索到一篇论文,里面详细介绍了云图的合成!

FY-4A多通道扫描辐射成像仪评价与图像合成
论文地址:http://www.doc88.com/p-8866426031707.html

python数字图像处理：图像数据类型及颜色空间转换:
https://www.cnblogs.com/denny402/p/5122328.html

绘制真彩图流程:

from netCDF4 import Dataset
import matplotlib.pyplot as plt
from skimage import io, data, img_as_float, img_as_ubyte, img_as_uint, img_as_int
import cv2hdf_data_path = "/Users/Downloads/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200805080000_20200805081459_4000M_V0001.HDF"nc_obj = Dataset(hdf_data_path)
type = nc_obj.variables.keys()
print(type)B = nc_obj.variables['NOMChannel01'][:]
G = nc_obj.variables['NOMChannel02'][:]
R = nc_obj.variables['NOMChannel03'][:]# 将数据类型转成int16
B = img_as_int(B)
G = img_as_int(G)
R = img_as_int(R)# opencv 将rgb三个通道融合
img3 = cv2.merge([R, G, B])img3 = exposure.adjust_log(img3, inv=True)#调整对比度
img3 = exposure.adjust_gamma(img3, 1.2)  # 图像调暗
T = 2
for i in range(len(img3)):for j in range(len(img3[0])):r = img3[i][j][0]g = img3[i][j][1]b = img3[i][j][2]img3[i][j] = (r * 0.6, g, b)if r / g > T:img3[i][j] = (g, r * 0.7, b)plt.imshow(img3, )
plt.axis('off')
plt.show()#

我这里只绘制出来的效果还是差点,最近比较忙有时间再解决吧!
上述问题,如果有解决了的同学,麻烦通知我一声

绘制中国地区的卫星云图

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import math
from numpy import deg2rad, rad2deg, arctan, arcsin, tan, sqrt, cos, sin
import numpy as np
from mpl_toolkits.basemap import Basemaphdf_data_path = "/Users/Downloads/FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200805080000_20200805081459_4000M_V0001.HDF"
ch_map = "/Users/map/中国地图/国界/bou1_4l"ea = 6378.137  # 地球的半长轴[km]
eb = 6356.7523  # 地球的短半轴[km]
h = 42164  # 地心到卫星质心的距离[km]
λD = deg2rad(104.7)  # 卫星星下点所在经度
# 列偏移
COFF = {"0500M": 10991.5,"1000M": 5495.5,"2000M": 2747.5,"4000M": 1373.5}
# 列比例因子
CFAC = {"0500M": 81865099,"1000M": 40932549,"2000M": 20466274,"4000M": 10233137}
LOFF = COFF  # 行偏移
LFAC = CFAC  # 行比例因子def latlon2linecolumn(lat, lon, resolution):"""经纬度转行列(lat, lon) → (line, column)resolution：文件名中的分辨率{'0500M', '1000M', '2000M', '4000M'}line, column"""# Step1.检查地理经纬度# Step2.将地理经纬度的角度表示转化为弧度表示lat = deg2rad(lat)lon = deg2rad(lon)# Step3.将地理经纬度转化成地心经纬度eb2_ea2 = eb ** 2 / ea ** 2λe = lonφe = arctan(eb2_ea2 * tan(lat))# Step4.求Recosφe = cos(φe)re = eb / sqrt(1 - (1 - eb2_ea2) * cosφe ** 2)# Step5.求r1,r2,r3λe_λD = λe - λDr1 = h - re * cosφe * cos(λe_λD)r2 = -re * cosφe * sin(λe_λD)r3 = re * sin(φe)# Step6.求rn,x,yrn = sqrt(r1 ** 2 + r2 ** 2 + r3 ** 2)x = rad2deg(arctan(-r2 / r1))y = rad2deg(arcsin(-r3 / rn))# Step7.求c,lcolumn = COFF[resolution] + x * 2 ** -16 * CFAC[resolution]line = LOFF[resolution] + y * 2 ** -16 * LFAC[resolution]return np.rint(line).astype(np.uint16), np.rint(column).astype(np.uint16)# 中国范围
x_min = 11
x_max = 54.75
y_min = 73.31
y_max = 135.91column = math.ceil((x_max - x_min) / 0.04)
row = math.ceil((y_max - y_min) / 0.04)
print(row, column)
ynew = np.linspace(y_min, y_max, row)  # 获取网格y
xnew = np.linspace(x_min, x_max, column)  # 获取网格x
xnew, ynew = np.meshgrid(xnew, ynew)  # 生成xy二维数组
data_grid = np.zeros((row, column))  # 声明一个二维数组keyword = "NOMChannel"nc_obj = Dataset(hdf_data_path)
type = nc_obj.variables.keys()
print(type)
print("--------------------------------------------")index = {}
r_data = {}
for k in type:if str(k).find(keyword) == 0:value = nc_obj.variables[k][:]for i in range(row):for j in range(column):lat = xnew[i][j]lon = ynew[i][j]fy_line = 0fy_column = 0if index.get((lat, lon)) == None:# 查找行列并记录下来下次循环使用fy_line, fy_column = latlon2linecolumn(lat, lon, "4000M")index[(lat, lon)] = fy_line, fy_columnelse:fy_line, fy_column = index.get((lat, lon))data_grid[i][j] = value[fy_line, fy_column]r_data[k] = data_gridimg = plt.figure()ax = img.add_subplot(111)m = Basemap(llcrnrlon=y_min, llcrnrlat=x_min, urcrnrlon=y_max, urcrnrlat=x_max)m.readshapefile(ch_map, 'states', drawbounds=True)p = plt.contourf(ynew, xnew, data_grid, cmap="gray", )plt.axis('off')plt.show()print("通道" + k + "绘制完成")

运行效果:
ok !!! 完事这样绘制出来的图片就可以叠加到地图上了。
代码看着简单,其实源文件里没让大家看,很多坑我这就不说了,希望大家直接避过坑,折腾我两个星期,不容易啊。

绘制出来的图片与官方对比

Python 解析风云四A卫星L1级别数据以及绘制卫星云图相关推荐

简单处理GPM数据和风云四号卫星数据
上一篇博客介绍了GPM数据和风云四号卫星数据的下载,本篇博客简介一些简单的处理和数据可视化操作. 一.GPM数据直接将下载好的文件直接拖入ArcMap中. 双击左侧图层下方的precipitatio ...
ENVI_IDL: 对风云四号卫星数据波段合成和线性拉伸并分别生成TIFF格式和JPEG格式
目录 1. 实验内容 2.编程 2.1 代码部分 2.2 运行结果 1. 实验内容 ---------------------------------------------------------- ...
[python]解析通达信盘后数据获取历史日线数据
平时我们在做离线的模型回溯测试时候,需要历史的k线数据. 可是通达信的日线数据如下: 日线数据在通达信的安装目录: vipdoc\sh\lday 下面本地的通达信是没有开放api和外部的 ...
python解析返回值类型为xml的数据接口
样参 <?xml version="1.0" encoding="UTF-8"?>-<root><state>success ...
Python解析罗永浩直播带货背后的数据秘密！
作为手机界最会说相声的罗永浩,已经正式加盟抖音,全身心投入直播行业了!按罗永浩的话说,是因为看了招商证券的调研报告,也为了偿还之前做手机留下来的债务,所以决定做电商直播.并且在4月1号晚间8点进行了直 ...
python解析tcp数据包-python解析获取发往本机的数据包并打印
1.[文件] tcp.py ~ 2KB 下载(69) # -*- coding: cp936 -*- import socket from struct import * from time impo ...
Python爬取‘跌妈不认’股票数据，绘制可视化图
前言本文的文字及图片来源于网络,仅供学习.交流使用,不具有任何商业用途,如有问题请及时联系我们以作处理. 以下文章来源于可以叫我才哥 ,作者才哥大家好,上次我们试着用vba在excel中绘制树状热 ...
python画多维散点图_python多维数据怎么绘制散点图
python matplotlib模块,是扩展的MATLAB的一个绘图工具库.他可以绘制各种图形,可是最近最的一个小程序,得到一些三维的数据点图,就学习了下python中的matplotlib模块,如 ...
python带货_Python解析罗永浩直播带货背后的数据秘密！
原标题:Python解析罗永浩直播带货背后的数据秘密! 作为手机界最会说相声的罗永浩,已经正式加盟抖音,全身心投入直播行业了!按罗永浩的话说,是因为看了招商证券的调研报告,也为了偿还之前做手机留下来的 ...
[系统安全] 四十一.APT系列(6)Python解析PE文件并获取时间戳判断来源区域
您可能之前看到过我写的类似文章,为什么还要重复撰写呢?只是想更好地帮助初学者了解病毒逆向分析和系统安全,更加成体系且不破坏之前的系列.因此,我重新开设了这个专栏,准备系统整理和深入学习系统安全.逆向分 ...

Python 解析风云四A卫星L1级别数据以及绘制卫星云图