Python 解析风云四A卫星L1级别数据以及绘制卫星云图 |
您所在的位置：网站首页 › 彩色绘画的地理图是什么 › Python 解析风云四A卫星L1级别数据以及绘制卫星云图 |
Python 解析风云四A卫星L1级别数据以及绘制卫星云图
|
绘制出来的卫星云图
全圆盘真彩图: 全圆盘单通道红外图:
数据格式说明:http://fy4.nsmc.org.cn/data/cn/data/realtime.html 数据下载地址:http://satellite.nsmc.org.cn/portalsite/Data/DataView.aspx?SatelliteType=1&SatelliteCode=FY4A# 本人使用的是4000M的全圆盘数据,下载数据需要申请账号 运行绘制出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 import math from numpy import deg2rad, rad2deg, arctan, arcsin, tan, sqrt, cos, sin import numpy as np from mpl_toolkits.basemap import Basemap hdf_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.求Re cosφe = cos(φe) re = eb / sqrt(1 - (1 - eb2_ea2) * cosφe ** 2) # Step5.求r1,r2,r3 λe_λD = λe - λD r1 = h - re * cosφe * cos(λe_λD) r2 = -re * cosφe * sin(λe_λD) r3 = re * sin(φe) # Step6.求rn,x,y rn = sqrt(r1 ** 2 + r2 ** 2 + r3 ** 2) x = rad2deg(arctan(-r2 / r1)) y = rad2deg(arcsin(-r3 / rn)) # Step7.求c,l column = 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.91 column = 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 = 0 fy_column = 0 if index.get((lat, lon)) == None: # 查找行列并记录下来下次循环使用 fy_line, fy_column = latlon2linecolumn(lat, lon, "4000M") index[(lat, lon)] = fy_line, fy_column else: fy_line, fy_column = index.get((lat, lon)) data_grid[i][j] = value[fy_line, fy_column] r_data[k] = data_grid img = 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 !!! 完事这样绘制出来的图片就可以叠加到地图上了。 代码看着简单,其实源文件里没让大家看,很多坑我这就不说了,希望大家直接避过坑,折腾我两个星期,不容易啊。
【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |