Python

这两个月可不好过啊，虽说进入了高三，当了学校的“老(xiao)大(di)哥(di)”，学业压力仍然不可避免的像一座山那样压在自己身上。

上面的话全部作废。因为我前面一个多个月实在是太浪了（还生病回家躺了半个月）。浪的结果就是月考（指9月两次考试）之后语文成绩直线下滑，数学重回比开学考还差水平（我不敢说下去了）。当然10月另说，毕竟进入了学习状态，除了地理英语“躺(zhuang)平(bi)”之外，10月月考基本上重回上学期的水平。

9月开始又开发了不少的东西，其中一个就是云图绘画脚本架构的更换。自从被别人（这个人是存在的，请在我以前的blog里的评论区找！）推荐使用Satpy模块之后，我开始了对Satpy模块的探究。用Satpy已经有一段时间了，然而当时那只是用于FY-4A云图的简单绘画当中，并未深入探索用法。直到我忍不住百度并且看(ying)了(ken)大半个小时文档，还和别人一起研究如何使用Satpy之后，终于知道了它真正的用法。

接下来我们将会介绍Satpy模块的使用方法。

本次我们使用linux CentOS 7.6.1810 x86_64环境，GOES ABI-L1b FLDK数据/FY4A AGRI FLDK数据/HIMAWARI-8 HSD TARGET AREA数据进行讲解

话不多说，先放一段关于satpy介绍的原文（请自行翻译）：

Satpy is a python library for reading, manipulating, and writing data from remote-sensing earth-observing meteorological satellite instruments. Satpy provides users with readers that convert geophysical parameters from various file formats to the common Xarray DataArray and Dataset classes for easier interoperability with other scientific python libraries. Satpy also provides interfaces for creating RGB (Red/Green/Blue) images and other composite types by combining data from multiple instrument bands or products. Various atmospheric corrections and visual enhancements are provided for improving the usefulness and quality of output images. Output data can be written to multiple output file formats such as PNG, GeoTIFF, and CF standard NetCDF files. Satpy also allows users to resample data to geographic projected grids (areas). Satpy is maintained by the open source Pytroll group.

其安装方式也在这个文档里面说得很清楚了：第一个是用conda install，第二个是Python pip。

'''Install methods from https://satpy.readthedocs.io/en/latest/install.html'''
# Conda-based Installation
$ conda config --add channels conda-forge
$ conda create -c conda-forge -n my_satpy_env python satpy
$ conda activate my_satpy_env
## In an existing environment
$ conda install -c conda-forge satpy

# Pip-based Installation
$ pip install satpy
$ pip install "satpy[viirs_sdr]"
$ pip install "satpy[all]"

# Ubuntu System Python Installation
$ sudo apt-get install python-pip python-gdal
$ sudo pip install virtualenv
$ virtualenv /path/to/pytroll-env
$ source /path/to/pytroll-env/bin/activate
$ pip install satpy

安装完之后，在终端上输入python，并输入import satpy，按下回车。如果没有报错，即为安装成功。

首先我们来一段正常人都看得懂（Who knows?）的Python代码（在这里，我们忽略了对satpy模块的配置。如果想要知道如何配置，请进入上文代码中的网址）：

import os, glob
from satpy import Scene

# First loading GOES-16 ABI-L1b data, and then create a Satpy Scene Object
filenames = glob.glob('/OR_ABI-L1b-RadF-M6C13_G16_s*_e*_c*.nc')
scn = Scene(filenames, reader='abi_l1b')

# Load C13 channel (Infrared channel)
channel = "C13"
scn.load([channel])

其中reader的值一定要对应你所读取的文件，如GOES的就用 abi_l1b 。以下是Satpy可用的reader列表：

正常人看其他地方的文档，会这么做：

# 保存到文件
# scn.save_dataset(channel, filename='{sensor}_{name}.png')

# 在notebook中显示
scn.show(channel)

虽然是画出来了，但是是绘画了全圆盘的数据，并且搞不好可能保存出来的图和你想象的不太一样（比如云顶黑的）。

当然如果你确实很懒，可以考虑直接用PIL模块甚至其它软件啥的直接反色。

至于往satpy做出来的图象里套入色阶，目前我还没从文档里看到对于色阶套入的教程，故不作进一步说明。（其实可以用matplotlib强行套一个，但是非常麻烦，本人没能完全做到）

那么如何截取区域云图呢？这里satpy/其他blog站点给出了案例（FY4A）：

# METHOD 01: Use pyresample to crop data
from pyresample import get_area_def

area_id = 'lekima'

x_size = 549
y_size = 499
area_extent = (-1098006.560556, -967317.140452, \
        1098006.560556, 1026777.426728)
projection = '+proj=laea +lat_0=19.0 +lon_0=128.0 +ellps=WGS84'
description = "Typhoon Lekima"
proj_id = 'laea_128.0_19.0'

areadef = get_area_def(area_id, description, proj_id, projection,x_size, y_size, area_extent)

# METHOD 02: Use coord2area_def.py to get area_def, and then crop data
$ python coord2area_def.py lekima_4km laea 10 28 118 138 4
lekima_4km:
  description: lekima_4km
  projection:
    proj: laea
    ellps: WGS84
    lat_0: 19.0
    lon_0: 128.0
  shape:
    height: 499
    width: 549
  area_extent:
    lower_left_xy: [-1098006.560556, -967317.140452]
    upper_right_xy: [1098006.560556, 1026777.426728]

# Then copy these text and save it, writing Python codes
import os, glob
from satpy import Scene

#  First loading FY-4A AGRI data, and then create a Satpy Scene Object
filenames = glob.glob('/home/FY4A-_AGRI*1000M_V0001.HDF')
scn = Scene(filenames, reader='agri_l1')

# load True Color composite channel
composite = 'true_color'
scn.load([composite])
scn.show(composite)


# if you were using METHOD 02, please write this code:
os.environ['PPP_CONFIG_DIR'] = '/home/satpy_config/'
lekima_scene = scn.resample('lekima_4km')

# if you were using METHOD 01, please write this code:
# lekima_scene = scn.resample(areadef)

lekima_scene.show(composite)
# lekima_scene.save_dataset(composite, filename='{sensor}_{name}_resampled.png')

然后出现以下图片的时候就说明成功啦：

使用pyresample进行satpy数据的截取只是一个非常简单的步骤。接下来我们将会讲解进阶方法（satpy+[basemap/matplotlib+cartopy]）。

这里以satpy+matplotlib+cartopy为例（博主无法安装basemap请谅解，如有配套basemap的方法可以在评论区发表）。

1.沿用上面的代码，照例载入数据至satpy中（数据采用HIMAWARI-8 TARGET AREA HSD）：

import os, glob
from satpy import Scene

# First loading GOES-16 ABI-L1b data, and then create a Satpy Scene Object
filenames = glob.glob('/home/HS_H08_*_B13_R30*')
scn = Scene(filenames, reader='abi_l1b')

# Load B13 channel (Infrared channel)
channel = "B13"
scn.load([channel])

2.获取加载数据的信息和经纬度数组等数据：

这里我就不得不说一下接下来要讲的截/获取数据的方法了。

在satpy里，当加载了数据之后，其数据内容通过reader被全部读取出来，保存至satpy.scene.Scene类之中。

而 satpy.scene.Scene类呢，（我自己也说不清是个啥东西）类似于一个xarray.DataArray对象，当然只在print或使用scn.to_xarray_dataset()时才会返回 xarray.DataArray 数据。读取数据之后的对象可以直接使用scn.value方法获取全数据（指对应通道的处理过的Brightness Temperature / Reflectance数据）。

# Get sensor finished scaning time
timestr = scn.end_time.strftime('%Y/%m/%d %H:%M:%SZ')

area_def = scn[channel].attrs['area']
lons, lats = area_def.get_lonlats()
data = scn[channel].values

# K unit to degC unit
data = data - 273.15

latmin, latmax, lonmin, lonmax = lats.mean() - 5, lats.mean() + 5, lons.mean() - 5, lons.mean() + 5

如果你担心全数据加载/绘图慢，satpy.scene.Scene类还有另外一种方法：scn.crop()

crop，说白了就是个裁取了数据范围之后的satpy.scene.Scene类，本质还是它。

下面是对satpy.scene.Scene.crop介绍的原文：

crop(area=None, ll_bbox=None, xy_bbox=None, dataset_ids=None)

Parameters

  • area (AreaDefinition) – Area to crop the current Scene to
  • ll_bbox (tuple, list) – 4-element tuple where values are in lon/lat degrees. Elements are (xmin, ymin, xmax, ymax) where X is longitude and Y is latitude.
  • xy_bbox (tuple, list) – Same as ll_bbox but elements are in projection units.
  • dataset_ids (iterable) – DataIDs to include in the returned Scene. Defaults to all datasets.

看完之后，再看看这段代码：

try:
    crop_scn = scn.crop(ll_bbox=(lon-0.5, lat-0.5, lon+0.5, lat+0.5))
except Exception:
    print("Satpy Error: Crop failed.")
    sys.exit(0)
else:
    area_def = crop_scn[channel].attrs['area']
    lons, lats = area_def.get_lonlats()
    data = crop_scn[channel].values
    data = data - 273.15
    # fix lower than 0 data to crop correctly
    lons[lons < 0] += 360

怎么样，看懂了么？可能会有人提问：为什么要用try-except Exception-else逻辑呢？为什么要调整经度数据呢？

前者是这样的：如果crop的区域超出了数据所在的范围，satpy会报错抛出个NotImplementedError，如果你想直接查看详细报错内容，可以去掉这个逻辑。

后者是这样的：当截取的经度超过180经度线时，即lon > 180，会自动把数据转为小于0的数值（即lon -= 360），当matplotlib绘图的时候，会根据经纬度数据进行对应位置的绘图，当有一部分大于0而一部分小于0的时候，会出现如下情况：

是不是感觉很莫名奇妙？但它确实如此。出现这种bug的原理可以看这个文档。

3.然后用matplotlib+cartopy进行绘图：

import matplotlib.pyplot as plt
import matplotlib.patches as mpathes
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

figsize = self.calc_figsize(georange)
dpi = 1200 / DEFAULT_WIDTH
background_color = "w"
f, ax = plt.subplots(figsize=figsize, subplot_kw=dict(projection=ccrs.PlateCarree(central_longitude=180)))
ax.set_extent([lonmin, lonmax, latmin, latmax], crs=ccrs.PlateCarree())
f.set_facecolor(background_color)
# close axis
plt.axis("off")
ax.add_feature(
    cfeature.COASTLINE.with_scale("10m"),
    facecolor="None",
    edgecolor="white",
    lw=0.5,
)
ax.add_feature(
    cfeature.LAKES.with_scale("10m"),
    facecolor="None",
    edgecolor="white",
    lw=0.25,
)
ax.add_feature(
    cfeature.RIVERS.with_scale("10m"),
    facecolor="None",
    edgecolor="white",
    linestyle='-',
    lw=0.25,
)
pcolor_kw = dict(cmap="gray_r", vmin=-80, vmax=40)
zorder = -1
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
gl = ax.gridlines(
    crs=ccrs.PlateCarree(),
    draw_labels=True,
    linewidth=0.6,
    linestyle=':',
    color='w',
    zorder=zorder
)
gl.rotate_labels = False
gl.xlabels_top = False
gl.xlabels_bottom = True
gl.ylabels_right = True
gl.ylabels_left = True
gl.xlabel_style = {'size': 4, 'color': 'k'}
gl.ylabel_style = {'size': 4, 'color': 'k'}
pcm = plt.pcolormesh(lons, lats, data, transform=ccrs.PlateCarree(), **pcolor_kw)
text = f"HIMAWARI-8 BAND13 TARGET AREA\nTime: {timestr}"
ax.set_title(text, loc='left', fontsize=6, color='k')
_ptext = f"{time}_{channel}"
plt.savefig(
    _ptext,
    dpi=dpi,
    bbox_inches="tight",
    pad_inches=0.03,
)
plt.close("all")

最后查看图片，出现类似于这样的图片就代表成功啦！

2021年第17号台风“狮子山”(Lionrock) 云图

更多关于Satpy的使用方法请浏览此文档，我们也可能在后续更新一些使用方法，敬请期待！

2 个评论

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注