风云卫星AWX格式读取与可视化
版本：Python3.7
数据：风云卫星AWX格式
数据来自：github
前言
近期需要读取awx格式数据，气象家园有人提到axw有对应的库，故测试一下awx的示例脚本
并作些简单美化
读取文件、访问数据、经纬度切片、转存netCDF文件
from awx import Awx
pathfile = '/home/mw/input/awx3540/ANI_IR2_R01_20230217_0800_FY2G.AWX'
ds = Awx(pathfile)
# 打印文件头信息
print(ds)
# 获取xarray.DataArray格式的观测数据
print(ds.values)
# 裁剪到指定的经纬度范围
print(ds.sel(lat=slice(20, 40), lon=slice(100, 130)))
# 保存数据为netCDF格式文件
ds.values.to_netcdf('ANI_VIS_R01_20230217_1000_FY2G.nc')
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1
AwxHead(filename_Sat96='ESLF170A.AWX', int_order=0, head1_length=40, head2_length=2112, filled_length=248, record_length=1200, head_record=3, data_record=1200, product_kind=1, zip_method=0, form_string='SAT2004', data_quality=0)
AwxGeosImageHead(satellite_name='FY2G', year=2023, month=2, day=17, hour=0, minute=0, channel=3, projection=1, width=1200, height=1200, ul_row_id=0, ul_pixel_id=0, sampling_rate=1, ul_lat=6206, lr_lat=659, ul_lon=7732, lr_lon=14870, clat=3500, clon=10000, std_lat1_or_lon=3000, std_lat2=6000, reso_h=500, reso_v=500, mark_geogrid_overlap=0, value_geogrid_overlap=255, palette_data_length=0, calibration_data_length=2048, position_data_length=0, reserved=0)
<xarray.DataArray (time: 1, y: 1200, x: 1200)>
array([[,
,
,
...,
,
,
]], dtype=uint8)
Coordinates:
* time (time) datetime64 2023-02-17
lat (y, x) float64 53.86 53.89 53.92 53.95 ... 5.944 5.933 5.922 5.912
lon (y, x) float64 50.12 50.18 50.25 50.31 ... 122.8 122.9 122.9 123.0
* x (x) float64 -3e+06 -2.995e+06 -2.99e+06 ... 2.995e+06 3e+06
* y (y) float64 3e+06 2.995e+06 2.99e+06 ... -2.995e+06 -3e+06
Attributes: (12/42)
filename_Sat96: ESLF170A.AWX
int_order: 0
head1_length: 40
head2_length: 2112
filled_length: 248
record_length: 1200
... ...
value_geogrid_overlap: 255
palette_data_length: 0
calibration_data_length:2048
position_data_length: 0
reserved: 0
proj4: +proj=lcc +lon_0=100.0 +lat_0=35.0 +lat_1=30.0 ...
<xarray.DataArray (time: 1, y: 588, x: 600)>
array([[[ nan,nan,nan, ...,nan,nan,nan],
[ nan,nan,nan, ...,nan,nan,nan],
[ nan,nan,nan, ...,nan,nan,nan],
...,
,
,
]])
Coordinates:
* time (time) datetime64 2023-02-17
lat (y, x) float64 46.63 46.63 46.63 46.63 ... 15.9 15.89 15.87 15.86
lon (y, x) float64 100.0 100.1 100.2 100.2 ... 126.0 126.1 126.1 126.1
* x (x) float64 2.502e+03 7.506e+03 1.251e+04 ... 2.995e+06 3e+06
* y (y) float64 1.254e+06 1.249e+06 1.244e+06 ... -1.679e+06 -1.684e+06
Attributes: (12/42)
filename_Sat96: ESLF170A.AWX
int_order: 0
head1_length: 40
head2_length: 2112
filled_length: 248
record_length: 1200
... ...
value_geogrid_overlap: 255
palette_data_length: 0
calibration_data_length:2048
position_data_length: 0
reserved: 0
proj4: +proj=lcc +lon_0=100.0 +lat_0=35.0 +lat_1=30.0 ...
利用Matplotlib绘制无投影的云图
import matplotlib.pyplot as plt
from awx import Awx
fpath = '/home/mw/input/awx3540/ANI_IR2_R01_20230217_0800_FY2G.AWX'
ds = Awx(pathfile=fpath)
print(ds)
dar = ds.values.squeeze()
plt.pcolormesh(dar.lon, dar.lat, dar, cmap='Greys_r')
plt.show()
AwxHead(filename_Sat96='ESLF170A.AWX', int_order=0, head1_length=40, head2_length=2112, filled_length=248, record_length=1200, head_record=3, data_record=1200, product_kind=1, zip_method=0, form_string='SAT2004', data_quality=0)
AwxGeosImageHead(satellite_name='FY2G', year=2023, month=2, day=17, hour=0, minute=0, channel=3, projection=1, width=1200, height=1200, ul_row_id=0, ul_pixel_id=0, sampling_rate=1, ul_lat=6206, lr_lat=659, ul_lon=7732, lr_lon=14870, clat=3500, clon=10000, std_lat1_or_lon=3000, std_lat2=6000, reso_h=500, reso_v=500, mark_geogrid_overlap=0, value_geogrid_overlap=255, palette_data_length=0, calibration_data_length=2048, position_data_length=0, reserved=0)
/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:8: UserWarning: The input coordinates to pcolormesh are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolormesh.
In :
fpath = '/home/mw/input/awx3540/FY2E_CTA_MLT_OTG_20170126_0130.AWX'
ds = Awx(pathfile=fpath)
print(ds)
dar = ds.values.squeeze()
plt.pcolormesh(dar.lon, dar.lat, dar, cmap='Greys_r')
plt.show()
AwxHead(filename_Sat96='DCZJ2613.AWX', int_order=0, head1_length=40, head2_length=80, filled_length=1081, record_length=1201, head_record=2, data_record=1201, product_kind=3, zip_method=0, form_string='SAT2004', data_quality=0)
AwxGridHead(satellite_name='FY2E', grid_element=20, grid_byte=1, grid_base=0, grid_scale=100, time_scale=0, start_year=2017, start_month=1, start_day=26, start_hour=1, start_minute=30, end_year=2017, end_month=1, end_day=26, end_hour=1, end_minute=55, ul_lat=6000, ul_lon=2700, lr_lat=-6000, lr_lon=14700, grid_unit=0, reso_h=10, reso_v=10, size_h=1201, size_v=1201, has_land=0, land=0, has_cloud=0, cloud=0, has_water=0, water=0, has_ice=0, ice=0, has_quality=1, quality_up=100, quality_down=0, reserve=0)
以数据指定的投影方式绘制云图
import os
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from awx import Awx
# fpath = r'./data/ANI_VIS_R02_20230217_1000_FY2G.AWX'# Mercator
fpath = '/home/mw/input/awx3540/ANI_IR2_R01_20230217_0800_FY2G.AWX'# lambert
ds = Awx(pathfile=fpath)
print(ds)
dar = ds.values.squeeze()
plt.figure(figsize=(8, 8))
if dar.projection == 1:
proj = ccrs.LambertConformal(central_longitude=dar.clon / 100,
central_latitude=dar.clat / 100,
standard_parallels=(dar.std_lat1_or_lon / 100.,
dar.std_lat2 / 100.))
extent =
elif dar.projection == 2:
proj = ccrs.Mercator(central_longitude=dar.clon / 100,
latitude_true_scale=dar.std_lat1_or_lon / 100.)
extent =
elif dar.projection == 4:
proj = ccrs.PlateCarree(central_longitude=dar.clon / 100.)
extent =
else:
raise NotImplementedError()
ax = plt.axes(projection=proj)
ax.set_extent(extent, crs=proj)
ax.coastlines(resolution='110m')
ax.gridlines(draw_labels=True)
ax.pcolormesh(dar.x, dar.y, dar, cmap='Greys_r')
plt.savefig(os.path.splitext(os.path.basename(fpath)) + '.png', dpi=300, bbox_inches='tight')
plt.show()
兰伯特投影经纬度优化版
用了WRF domain 绘制改进的经纬度函数，这里就不复制粘贴了。
import osdd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from awx import Awx
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import shapely.geometry as sgeom
import numpy as np
import matplotlib.transforms as mtransforms
import cartopy.crs as ccrs
from cartopy.io.shapereader import BasicReader
import cmaps
from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
from mpl_toolkits.axes_grid1.inset_locator import TransformedBbox, BboxPatch, BboxConnector
import shapely.geometry as sgeom
from copy import copy
import concurrent.futures
fpath = '/home/mw/input/awx3540/ANI_IR2_R01_20230217_0800_FY2G.AWX'# lambert
provinces = BasicReader('/home/mw/input/data5246/中国地图/China_provinces/China_provinces.shp')
countries = BasicReader('/home/mw/input/data5246/中国地图/world_countries/world_countries.shp')
ds = Awx(pathfile=fpath)
print(ds)
dar = ds.values.squeeze()
if dar.projection == 1:
proj = ccrs.LambertConformal(central_longitude=dar.clon / 100,
central_latitude=dar.clat / 100,
standard_parallels=(dar.std_lat1_or_lon / 100.,
dar.std_lat2 / 100.))
extent =
elif dar.projection == 2:
proj = ccrs.Mercator(central_longitude=dar.clon / 100,
latitude_true_scale=dar.std_lat1_or_lon / 100.)
extent =
elif dar.projection == 4:
proj = ccrs.PlateCarree(central_longitude=dar.clon / 100.)
extent =
else:
raise NotImplementedError()
fig = plt.figure(figsize=(10, 8),dpi=100)
ax = plt.axes(projection=proj)
ax.set_extent(extent, crs=proj)
# 添加经纬度网格和刻度
# *must* call draw in order to get the axis boundary used to add ticks:
fig.canvas.draw()
# Define gridline locations and draw the lines using cartopy's built-in gridliner:
xticks=
yticks=
ax.gridlines(xlocs=xticks, ylocs=yticks)
font3={'family':'SimHei','size':14,'color':'k'}
ax.gridlines(xlocs=xticks, ylocs=yticks, draw_labels=False, linewidth=0.8, color='k', alpha=0.6, linestyle='--')
# Label the end-points of the gridlines using the custom tick makers:
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER)
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)
ax.pcolormesh(dar.x, dar.y, dar, cmap='Greys_r')
ax.add_geometries(provinces.geometries(), linewidth=.5, edgecolor='white', crs=ccrs.PlateCarree(),
facecolor='none')
ax.add_geometries(countries.geometries(), linewidth=.5, edgecolor='black', crs=ccrs.PlateCarree(),
facecolor='none')
plt.show()
简单易懂，上手无门槛
文章来源于微信公众号：气ython风雨
有没有针对初学者的更简单的解释？
页:
[1]