风云卫星AWX格式读取与可视化

用户11172986

发布于 2024-06-20 16:55:54

700

发布于 2024-06-20 16:55:54

文章被收录于专栏：气python风雨气python风雨

版本：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([[[202, 202, 203, ..., 185, 185, 185],
        [201, 202, 202, ..., 185, 184, 184],
        [201, 201, 201, ..., 184, 184, 184],
        ...,
        [109, 109, 110, ..., 128, 127, 126],
        [109, 109, 109, ..., 127, 126, 126],
        [109, 109, 109, ..., 125, 125, 125]]], dtype=uint8)
Coordinates:
  * time     (time) datetime64[ns] 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],
        ...,
        [190., 192., 192., ...,  nan,  nan,  nan],
        [189., 191., 191., ...,  nan,  nan,  nan],
        [187., 189., 190., ...,  nan,  nan,  nan]]])
Coordinates:
  * time     (time) datetime64[ns] 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 [4]:

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 = [dar.x.min(), dar.x.max(), dar.y.min(), dar.y.max()]
elif dar.projection == 2:
    proj = ccrs.Mercator(central_longitude=dar.clon / 100,
                         latitude_true_scale=dar.std_lat1_or_lon / 100.)
    extent = [dar.x.min(), dar.x.max(), dar.y.min(), dar.y.max()]
elif dar.projection == 4:
    proj = ccrs.PlateCarree(central_longitude=dar.clon / 100.)
    extent = [dar.lon.min(), dar.lon.max(), dar.lat.min(), dar.lat.max()]
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))[0] + '.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 = [dar.x.min(), dar.x.max(), dar.y.min(), dar.y.max()]
elif dar.projection == 2:
    proj = ccrs.Mercator(central_longitude=dar.clon / 100,
                         latitude_true_scale=dar.std_lat1_or_lon / 100.)
    extent = [dar.x.min(), dar.x.max(), dar.y.min(), dar.y.max()]
elif dar.projection == 4:
    proj = ccrs.PlateCarree(central_longitude=dar.clon / 100.)
    extent = [dar.lon.min(), dar.lon.max(), dar.lat.min(), dar.lat.max()]
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=[80,85,90,95,100,105,110,115,120]
yticks=[10,15,20,25,30,35,40,45,50,55]
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()