目录（暂定）
1. 给初学者的建议
2. 卫星数据读取——以MODIS为例
a. MODIS数据下载
b. Python环境
c. HDF文件读取，MODIS波段简介
d. 辐射定标与亮温反演
e. 经纬度插值
3. 红外云图渲染
a. 色阶
b. Matplotlib出图
4. 伪VIS云图渲染
a. 伪VIS基础算法实现
b. 海温反演修正与高云修正

卫星数据读取——以MODIS为例
a. MODIS数据下载
MODIS数据可以从以下两个网站获取：
1. Earthdata search https://search.earthdata.nasa.gov/search
2. LAADS DAAC https://ladsweb.modaps.eosdis.nasa.gov/search/order
追风时画实况云图我会偏向于使用Earthdata，画历史云图有时候会使用LAADS，云图预览方便些，这里我们主要介绍第一个网站。

这个网站目前提供了9331种不同的数据集，其中绝大部分是我们所不感兴趣的。左侧筛选栏中选择Instruments-MODIS，Processing Levels-1B，就可以筛选出MODIS L1B的数据，一共21个，当然这里还包含一些我们不感兴趣的数据（比如5km重采样）：

MODIS数据集有各自的代号，例如MYD021KM是Aqua上的1km分辨率L1B数据，其中Y/O代表卫星Aqua/Terra，1KM/HKM/QKM代表1km/0.5km/0.25km分辨率。（H=Half, Q=Quarter）
本贴主要使用的是1km分辨率的MYD021KM/MOD021KM数据，文件中包含MODIS所有波段的数据；由于不是所有波段都有0.5km与0.25km的分辨率，HKM和QKM仅包含部分波段的数据，绘制VIS和真彩色时可能用到。
网站中部分数据是NRT（Near Real Time）版本，更新相对更快（一般在卫星通过后1~3小时不等出数据）但只有近期存档，数据处理方式和正常版本一致，适合实时追风时使用。
选择数据集后可以点击左上角的小图标筛选指定时间和矩形区域的数据：

（大部分卫星数据网站都会提供类似筛选数据的功能，也有少数几个得自己对着时间翻找的，比如ncei的微波存档）
近期Earthdata增加了上图右侧的云图预览，但个人测试下来这个功能会将每颗卫星白天和夜晚的两次扫描混淆，暂时不建议使用。左侧数据自带的预览是正常可信的，可以用于确认目标TC是否扫正。
这里我们选中了海燕11.071345Z的Terra数据并下载。境内访问Earthdata有可能比较卡，但这个问题需要各位自己解决。下载过程如果非常慢可以尝试用IDM等工具加速。

前面提到过绘制云图有时候也需要参考一些文献。MODIS L1B数据的文献可以在这里获取：https://mcst.gsfc.nasa.gov/l1b-documents
官方给的文献一般比较长（比如这里User Guide有六十多页），不需要全部读完，靠目录或者关键词搜索找到自己想要的信息就行了。后续HDF文件读取与辐射定标等步骤均可以在User Guide里面找到对应的说明，教程里将会以更简单的方式呈现这些步骤，但如果想要了解“为什么这样处理数据”，去翻一翻UserGuide也是不错的选择。

b. Python环境
Python环境分为开发环境和第三方库环境。合适的开发环境可以让代码页面更清晰美观，区分各种类型的变量，检查部分显而易见的bug，并提供一些（暂时不一定用得上的）辅助功能。
Python本体与开发环境的安装大家可以参考bala近期写的Python处理气象数据教程的第一、二章（https://www.tyboard.net/forum.php?mod=viewthread&tid=316），我个人开发环境选的是Pycharm，总体而言用的还算顺手，就是内存占的比较大，具体选哪个开发环境大家可以依自己需求决定。
成熟、易于使用的第三方库是Python的一大优势。绘制云图时我们无需将最底层的代码全部自己码完（很显然我也不知道这怎么写），只需要学会调用第三方库中已经编写好的模块就行。
完成最基础的MODIS云图渲染只需要调用四个第三方库：pyhdf，numpy，scipy，matplotlib，分别负责hdf读取，数学操作，数据插值（当然scipy能做的事远不止插值），与图片渲染。如果有绘制海岸线等需求会用到cartopy，伪VIS使用陆地底图还会用到opencv等。
第三方库的安装同样可以参考bala的教程的第四章，（https://www.tyboard.net/forum.php?mod=viewthread&tid=316）写的很详细，我自己写不了这么周全就不作重复了。使用pip安装时请预先确认第三方库在pip中对应的名称，有时会和调用代码时用的名称不一致。
大型的第三方库一般都有网站对其功能作详细的解释。以matplotlib为例，官网上（https://matplotlib.org/）自带了一些代码示例。不过以个人体验来说，更多时候我都是直接在浏览器上搜函数名称，跳转到matplotlib.org的reference页面去查看函数具体的输入输出，若理解有困难再去参考页面底部链接的代码示例。


上图是pcolormesh函数的reference页面，这也是我们将亮温等二维数组绘制成投影过的云图最主要的函数。

c. HDF文件读取，MODIS波段简介
卫星云图数据存储的格式相当多样，其中不少格式需要较复杂的脚本才能读取（比如Himawari-8/9最广泛使用的DAT格式）。相对容易读取的数据存储方式包括HDF4，HDF5，NetCDF等。MODIS的L1B数据正是用HDF4格式存储的。HDF4有多种读取方式，这里我们选择使用第三方库pyhdf。HDF5和NetCDF我使用的分别是h5py和netCDF4，在此不作详细介绍。
pyhdf的介绍：https://fhs.github.io/pyhdf/modules/SD.html，下图列举了部分会用到的功能

from pyhdf.SD import SD, SDC
import numpy as np
首先导入pyhdf与numpy，此处pyhdf只导入了其中的一部分模块，numpy则是导入时命名为了np，后续反复调用的时候少打几个字母。
file = SD(r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf', SDC.READ)
在只读模式（SDC.READ）下打开文件，文件路径请自行替换，这里是用了Patricia10.23 0515z的Terra，此处路径前的r是为了直接拷贝过来的绝对路径能够正常读取。
我们先使用datasets()看看文件中放了哪些数据，这一步相当于列出数据集的“目录”，还没有真正将数据取出来：
data = file.datasets()
print(data)
print(data.keys())
直接print(data)会返回每个数据集的名称与一些附带信息，加上keys()能仅返回数据集名称，以下是print(data.keys())的结果：
dict_keys(['Latitude', 'Longitude', 'EV_1KM_RefSB', 'EV_1KM_RefSB_Uncert_Indexes', 'EV_1KM_Emissive', 'EV_1KM_Emissive_Uncert_Indexes', 'EV_250_Aggr1km_RefSB', 'EV_250_Aggr1km_RefSB_Uncert_Indexes', 'EV_250_Aggr1km_RefSB_Samples_Used', 'EV_500_Aggr1km_RefSB', 'EV_500_Aggr1km_RefSB_Uncert_Indexes', 'EV_500_Aggr1km_RefSB_Samples_Used', 'Height', 'SensorZenith', 'SensorAzimuth', 'Range', 'SolarZenith', 'SolarAzimuth', 'gflags', 'EV_Band26', 'EV_Band26_Uncert_Indexes', 'Band_250M', 'Band_500M', 'Band_1KM_RefSB', 'Band_1KM_Emissive', 'Noise in Thermal Detectors', 'Change in relative responses of thermal detectors', 'DC Restore Change for Thermal Bands', 'DC Restore Change for Reflective 250m Bands', 'DC Restore Change for Reflective 500m Bands', 'DC Restore Change for Reflective 1km Bands'])
其中'Latitude'，'Longitude'是经纬度数据，'EV_1KM_RefSB'是反射波段数据，这里暂时用不到，'EV_1KM_Emissive'是热辐射波段数据，包含了我们所需要的各个红外波段。
使用select将'Latitude'，'Longitude'，'EV_1KM_Emissive'读取出来，并用get()读取数组，用np.shape()获取每个数组的形状（不过其实print(data)也会将数组形状作附带信息返回）：
lon = file.select('Longitude').get()
lat = file.select('Latitude').get()
Emissive = file.select('EV_1KM_Emissive').get()
print(lon, lat, Emissive)
print(np.shape(lon), np.shape(lat), np.shape(Emissive))
可以看到热辐射波段数组的形状是(16, 2030, 1354)，代表有16个不同波段堆叠在一起，每个波段的图像有2030×1354个像素，而数组的值都是0~65535的整数，显然不是我们最终需要的亮温，需经过辐射定标与亮温反演（见2.d）；经纬度数组的值是正确的，但数组形状是(406, 271)，像素比波段数据少，这是因为文件中存储的经纬度降过分辨率，需插值后才能使用（见2.e）。
'EV_1KM_Emissive'这个数据集并非只包含上面读取的数组，HDF文件可以为数据集附上attributes，其中往往包含辐射定标的关键信息，查看方式如下：
print(file.select('EV_1KM_Emissive').attributes())
print(file.select('EV_1KM_Emissive').attributes().keys())
以下是加上keys()，只获取attributes名称的结果：
dict_keys(['long_name', 'units', 'valid_range', '_FillValue', 'band_names', 'radiance_scales', 'radiance_offsets', 'radiance_units'])
其中'long_name'是数据集全称， 'band_names'是波段号，这两个attribute可以在看官方文献看晕掉的时候让你知道读的是什么数据；'radiance_scales'，'radiance_offsets'是辐射定标的两个系数（见2.d），'radiance_units'是辐射定标后得到的辐亮度的单位，有时候不同卫星数据会采用不同的辐亮度单位，反演的亮温如果偏差很大则需要检查该单位并调整黑体辐射公式的形式。
获取attributes中的数组同样是使用get()即可：
radiance_scales = file.select('EV_1KM_Emissive').attributes().get('radiance_scales')
radiance_offsets = file.select('EV_1KM_Emissive').attributes().get('radiance_offsets')
MODIS波段介绍：https://modis.gsfc.nasa.gov/about/specifications.php
先前查看热辐射波段attributes的时候16个波段号分别为20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36。细心的读者可能已经发现了，Band26不在其中。（不细心的作者因为这个在伪VIS教程篇里面把中波红外要用的Band20错写成Band21了，前两天刚勘误）Band26是短波红外反射波段，在2.d中读取指定波段数据的时候需要注意到这一点，不然可能会读错波段。
16个热辐射波段中，Band31用于渲染我们最熟知的IR云图，Band27则是标准水汽波段，Band28是响应层面偏低一些的水汽波段，反映中低层大气的暖心和水汽吸收情况。伪VIS所用到的波段则是Band31，Band32与Band20，波段选择的具体原因可以参考《伪VIS算法原理与教程 - 理论篇》。值得一提的是，Band20与Band21,22的波长非常相近，选择Band20的主要原因是Band20接收的波长宽一些，更容易接收到较冷云顶上微弱的辐射，云顶噪点相对而言少一些。

d. 辐射定标与亮温反演
辐射定标将卫星仪器测量的原始数据转化为辐亮度这个物理量，而亮温反演通过黑体辐射公式将辐亮度转化为亮温。
从卫星原始数据计算出辐亮度是个很复杂的过程，需要考虑卫星本身的电子元件问题等等，好在L1B数据已经将这些效应处理掉了，只需要我们做辐射定标的最后一步，把波段数据里存储的整数转化为辐亮度。MODIS User guide（https://mcst.gsfc.nasa.gov/sites/default/files/file_attachments/M1054E_PUG_2022_1005_V6.2.2_Terra_V6.2.3_Aqua.pdf）中的equation 5.8给出了辐射定标的公式：

其中SI（Scaled Integer）就是波段数据存储的整数，L则是辐亮度。
具体到python代码上：
radiance = (Emissive[10] - radiance_offsets[10]) * radiance_scales[10]
这里计算了Band31的辐亮度。波段数据和辐射定标系数均按照波段序号排列，跳过Band26这个反射波段，Band31是热辐射波段中第11个波段，而python索引从0开始算，因此索引[10]正对应Band31。
黑体辐射的概念在《伪VIS算法原理与教程》理论篇2.1章节有较为详细的介绍，这里只放由此而来的亮温计算公式：

其中T是开氏度单位的亮温，转换到摄氏度需要-273.15，ℎ是普朗克常数，𝑐是光速，𝑘是玻尔兹曼常数，𝜆是波长，ln()是取自然对数，B𝜆是辐亮度。对常数做一下简化：
带入计算可得在SI制单位下C1=0.0143877，C2=1.191043*10^-16，但红外波段的波长一般用微米表示，且查看attributes（见2.c）可知辐亮度也用到了微米作单位，因此调整单位后C1=1.43877*10^4，C2=1.191043*10^8。
带入Band31中心波长（11.030微米）即可将亮温计算出来，其中**代表指数运算，np.log取的就是自然对数：
C1 = 1.4387685*10**4
C2 = 1.19104356*10**8
bt31 = C1/11.030 / (np.log(1 + C2/11.030**5/radiance)) - 273.15
至此我们已经获得了MODIS Band31的准确摄氏度亮温。

贴吧移动端最近增加了一个隐藏/展开全文的特性，但看起来还不完善，我这里会出现段落重复，需要点击全文才能正常显示

e. 经纬度插值
先前提到过MODIS L1B文件中的经纬度数值降过分辨率，云图像素数为2030×1354，而经纬度像素数为406×271。为了正常投影云图，我们需要将经纬度数组插值到与云图像素一一对应。我们可以用 scipy.interpolate.RectBivariateSpline类 实现插值。
RectBivariateSpline类，提供样条插值功能：

导入scipy.interpolate
import scipy.interpolate
首先提取出经度的形状406×271
ny, nx = lon.shape
用np.arange创建两个数值为0~405，0~270，间隔为1的一维数组，这代表经纬度原始的矩形网格
x = np.arange(nx)
y = np.arange(ny)
再创建数值为0~405，0~270，间隔为271/1354，1/5的一维数组，数组长度为2030，1354，这代表插值后的矩形网格
newx = np.arange(0, nx, 271/1354)
newy = np.arange(0, ny, 1/5)
应用RectBivariateSpline类，首先明确原网格上每个点的经度数值，随后用刚生成的数组将其插值到2030×1354的网格上
splinelon = scipy.interpolate.RectBivariateSpline(y, x, lon)
newlon = splinelon(newy, newx)
此处newlon就是插值后的经度，要获取插值后的纬度重复最后一个步骤并用lat代替即可。有了一一对应的亮温、经度、纬度三个数组，我们就可以开始绘制投影后的红外云图了。

红外云图渲染
a. 色阶
相信各位风迷对色阶并不陌生，色阶将红外云图的每个像素按照其亮温上色。色阶颜色随亮温的变化可以分为渐变（如WV色阶0度以下的部分）和跳变（如IR-BD OW以下的部分），当然许多色阶两者兼有（如IR-RAMMB），黑白的IR和VIS云图本质上也可以看做使用了一类从黑到白的渐变色阶。

左侧是Patricia10.23 0515z Terra Band31未经投影的IR-Rammb图像，右侧是IR-Rammb，IR-BD，WV，IR-BW的色阶条。
色阶编写需要用到matplotlib.colors.LinearSegmentedColormap类（https://matplotlib.org/stable/ap ... mentedColormap.html），最简单的IR-BW色阶长这样：
bwdata={'green': [(0, 1, 1),
(1, 0, 0)],
'red': [(0, 1, 1),
(1, 0, 0)],
'blue': [(0, 1, 1),
(1, 0, 0)]}
卫星云图每个像素的颜色可以用RGB值表示，如果想要自己编写色阶，可以拿读取像素RGB值的软件找找感觉。LinearSegmentedColormap的色阶格式也有'green' 'red' 'blue'独立的三组，对应G，R，B的亮度值。IR-BW是黑白色阶，因此三组都是一样的，彩色色阶则会有所不同。
每一行色阶由三个数字组成（称为x, y0, y1），并指定某一亮温节点对应的该颜色（RGB其中之一）亮度。x代表该亮温节点的位置，当色阶渲染范围是-100到50度时，x为0代表-100度，为0.5代表-25度，以此类推。y0, y1代表该颜色亮度，其0~1的范围对应0~255的RGB值。为什么需要用两个数字代表同一个亮温节点上的亮度呢？这是因为每个亮温节点两侧的亮度可以不一样，出现颜色的跳变。y0代表更低亮温那一侧的亮度，y1代表更高亮温那一侧的亮度。在两个亮温节点的间隙中，LinearSegmentedColormap会根据上一行的y1与下一行的y0自动进行渐变，若这两个值相同则无渐变。对于上面的IR-BW色阶来说，就是在-100度 亮度为1 和50度 亮度为0的区间内进行渐变。

下面展示几个更复杂的色阶实例：
wvdata={'green': [(0, 0/255, 0/255),
(30/150, 0/255, 0/255),
(45.5/150, 128/255, 128/255),
(58.5/150, 255/255, 255/255),
(65/150, 255/255, 255/255),
(72.5/150, 255/255, 255/255),
(86/150, 128/255, 128/255),
(100/150, 20/255, 255/255),
(1, 1, 1)],
'red': [(0, 128/255, 128/255),
(30/150, 128/255, 128/255),
(45.5/150, 255/255, 255/255),
(58.5/150, 255/255, 255/255),
(65/150, 128/255, 128/255),
(72.5/150, 128/255, 128/255),
(86/150, 0/255, 0/255),
(100/150, 100/255, 255/255),
(1, 1, 1)],
'blue': [(0, 0/255, 0/255),
(30/150, 0/255, 0/255),
(45.5/150, 0/255, 0/255),
(58.5/150, 128/255, 128/255),
(65/150, 128/255, 128/255),
(72.5/150, 255/255, 255/255),
(86/150, 255/255, 255/255),
(100/150, 100/255, 255/255),
(1, 1, 1)]}
WV色阶应用于-100~50度的渲染范围，因此x0位置上的N/150就对应-100+N度。WV色阶以渐变为主，因此前几行色阶y0y1是相同的，代表没有颜色的跳变，上一行y1与下一行y0不同，代表在亮温节点之间进行渐变。注意在x0=100/150也就是0度的时候y0y1不同，代表此处出现了从紫色到白色的跳变，这一点可以通过本节第一张图确认。
bddata = {'blue': [(0.0, 0.0, 0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24, 0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0), (1.0, 0.0, 0.0)],
'green': [(0.0, 0.0, 0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24,0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0),
(1.0, 0.0, 0.0)],
'red': [(0.0, 0.0,0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24, 0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0),
(1.0, 0.0, 0.0)]}
BD色阶中有些小数没化为分数，不影响使用，但确实不太美观（
BD色阶同样应用于-100~50度的渲染范围。仔细观察前几行的y0，y1不难发现，同一行的y0，y1不同，代表亮温节点上会出现颜色的跳变，例如CDG在-81跳变到CMG；上一行y1与下一行y0相同，代表亮温节点之间不存在渐变，对应着从-81~-76区间的任意亮温都会被赋予CMG的颜色。

3.红外云图渲染a.色阶这节被系统删了，先看看能不能恢复吧，不然发个截图也不方便拷贝色阶代码不过去论坛看也行

b. Matplotlib出图
在之前的步骤中我们已经得到了云图每个像素的亮温，经纬度，以及用于上色的色阶，下一步就是将云图真正渲染出来。
matplotlib.pyplot.imshow是二维数组可视化最简单的方式之一，我们可以直接把Band31亮温这个数组画出来：
import matplotlib.pyplot as plt
plt.imshow(bt31)
plt.savefig('文件名')
以Patricia为例，运行代码会在.py的同目录下生成如下图片。扫描边缘有拉伸且方向有反转是未经过投影导致的。此处没有设置色阶，因此用的是matplotlib的默认蓝黄色阶。

通过cmap参数设置一下bd色阶
plt.imshow(bt31, cmap=LinearSegmentedColormap('1', bddata))

看着确实是BD色阶，但很明显Patricia没有这么多CDG。色阶的偏差是未设置色阶渲染范围导致的，这会使matplotlib自动取数组的最大最小值作为渲染范围。为解决这个问题我们加上vmin，vmax固定渲染范围：
plt.imshow(bt31, vmin=-100, vmax=50, cmap=LinearSegmentedColormap('1', bddata))

投影需要另一个函数pcolormesh完成，使用的参数和imshow差不多，只是多了（插值后的）经度，纬度两个数组。这里默认是等经纬度投影。
plt.pcolormesh(newlon, newlat, bt31, vmin=-100, vmax=50, cmap=LinearSegmentedColormap('1', bddata))

上述投影方式会在日界线附近经度从180跳到-180时失效，这里用了个相当凑合的解决方法，将所有负值的经度值+360，使其在日界线附近连续，但在本初子午线附近不连续：
lon[lon < 0] += 360
前面几张图片的大小显然不是我们想要的，完全没有发挥出MODIS应有的分辨率。图片大小可以由plt.figure的figsize参数设置，这里贴一下我比较常用的设置方法，不过只适用于等经纬度投影，肯定不是最好的方法但比较简单：
londiff = np.max(lon)-np.min(lon) #计算图像两端经纬度差
latdiff = np.max(lat)-np.min(lat)
plt.figure(figsize=(londiff, latdiff), dpi=150) #按照图像两端经纬度差设置像素数，每经纬度150像素，可按需调整
plt.subplots_adjust(top=1, bottom=0, left=0, right=1, hspace=0, wspace=0) #使图片充满画布
plt.axis('off') #关闭坐标轴
加上这些设置得到的图片如下，和MODIS存档贴里我和qc等人更新的没什么区别了：

应读者要求贴一个非常简短的demo，其中每一行在教程中都有比较详细的介绍。
demo可以直接运行以渲染MODIS Band31红外云图，不过暂不包含其他波段以及数据切片，读眼温等功能，这些稍后会单独讲一下。
主文件：
file = SD(r'E:\Data\MOD021KM.A2015296.0515.061.2017322234443.hdf', SDC.READ) #2.c HDF文件读取
data = file.datasets()
lon = file.select('Longitude').get()
lat = file.select('Latitude').get()
Emissive = file.select('EV_1KM_Emissive').get()
radiance_scales = file.select('EV_1KM_Emissive').attributes().get('radiance_scales')
radiance_offsets = file.select('EV_1KM_Emissive').attributes().get('radiance_offsets')
radiance = (Emissive[10] - radiance_offsets[10]) * radiance_scales[10] #2.d 辐射定标与亮温反演
C1 = 1.4387685*10**4
C2 = 1.19104356*10**8
bt31 = C1/11.030 / (np.log(1 + C2/11.030**5/radiance)) - 273.15
ny, nx = lon.shape #2.e 经纬度插值
x = np.arange(nx)
y = np.arange(ny)
newx = np.arange(0, nx, 271/1354)
newy = np.arange(0, ny, 1/5)
splinelon = scipy.interpolate.RectBivariateSpline(y, x, lon)
newlon = splinelon(newy, newx)
splinelat = scipy.interpolate.RectBivariateSpline(y, x, lat)
newlat = splinelat(newy, newx)
londiff = np.max(lon)-np.min(lon) #3.b Matplotlib出图
latdiff = np.max(lat)-np.min(lat)
plt.figure(figsize=(londiff, latdiff), dpi=150)
plt.subplots_adjust(top=1, bottom=0, left=0, right=1, hspace=0, wspace=0)
plt.axis('off')
plt.pcolormesh(newlon, newlat, bt31, vmin=-100, vmax=50, cmap=LinearSegmentedColormap('1', bddata))
plt.savefig('文件名')
色阶文件（cmap.py，放在主文件同目录下）：
bddata = {'blue': [(0.0, 0.0, 0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24, 0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0), (1.0, 0.0, 0.0)],
'green': [(0.0, 0.0, 0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24,0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0),
(1.0, 0.0, 0.0)],
'red': [(0.0, 0.0,0.3333333333333333),
(0.12666666666666668, 0.3333333333333333, 0.5294117647058824),
(0.16, 0.5294117647058824, 1.0),
(0.2, 1.0, 0.0),
(0.24, 0.0, 0.6274509803921569),
(0.30666666666666664, 0.6274509803921569, 0.43137254901960786),
(0.38666666666666666, 0.43137254901960786, 0.23529411764705882),
(0.46, 0.23529411764705882, 0.792156862745098),
(0.7266666666666667, 0.42745098039215684, 1.0),
(0.8533333333333334, 0.0, 0.0),
(1.0, 0.0, 0.0)]}
rammbdata={'red': [(0.0, 0.0, 0.3137254901960784),
(0.06666666666666667, 1.0, 0.3137254901960784),
(0.13333333333333333, 1.0, 1.0),
(0.2, 0.39215686274509803, 0.0),
(0.26666666666666666, 0.0, 0.0),
(0.3333333333333333, 0.0, 0.3333333333333333),
(0.4666666666666667, 0.7058823529411765, 1.0),
(0.9, 0.0, 0.0), (1.0,0.0, 0.0)],
'blue': [(0.0, 0.0, 0.3137254901960784),
(0.06666666666666667, 1.0,0.3137254901960784),
(0.13333333333333333, 0.0, 0.0),
(0.2, 0.0, 0.0),
(0.26666666666666666, 0.0, 1.0),
(0.3333333333333333, 0.39215686274509803, 0.3333333333333333),
(0.4666666666666667, 1.0, 1.0),
(0.9, 0.0, 0.0),
(1.0, 0.0, 0.0)],
'green': [(0.0, 0.0, 0.3137254901960784),
(0.06666666666666667, 1.0, 0.3137254901960784),
(0.13333333333333333, 1.0, 0.0),
(0.2, 0.0, 1.0),
(0.26666666666666666, 0.39215686274509803, 0.0),
(0.3333333333333333, 0.0, 0.3333333333333333),
(0.4666666666666667, 1.0, 1.0),
(0.9, 0.0, 0.0),
(1.0, 0.0, 0.0)]}
wvdata={'green': [(0, 0/255, 0/255),
(30/150, 0/255, 0/255),
(45.5/150, 128/255, 128/255),
(58.5/150, 255/255, 255/255),
(65/150, 255/255, 255/255),
(72.5/150, 255/255, 255/255),
(86/150, 128/255, 128/255),
(100/150, 20/255, 255/255),
(1, 1, 1)],
'red': [(0, 128/255, 128/255),
(30/150, 128/255, 128/255),
(45.5/150, 255/255, 255/255),
(58.5/150, 255/255, 255/255),
(65/150, 128/255, 128/255),
(72.5/150, 128/255, 128/255),
(86/150, 0/255, 0/255),
(100/150, 100/255, 255/255),
(1, 1, 1)],
'blue': [(0, 0/255, 0/255),
(30/150, 0/255, 0/255),
(45.5/150, 0/255, 0/255),
(58.5/150, 128/255, 128/255),
(65/150, 128/255, 128/255),
(72.5/150, 255/255, 255/255),
(86/150, 255/255, 255/255),
(100/150, 100/255, 255/255),
(1, 1, 1)]}
bwdata={'green': [(0, 1, 1),
(1, 0, 0)],
'red': [(0, 1, 1),
(1, 0, 0)],
'blue': [(0, 1, 1),
(1, 0, 0)]}

简单提几个绘制红外云图可能会用到的杂项功能，暂时不作特别详细的展开，以后可能会在本楼慢慢补充/勘误：
1. 云图裁剪
MODIS的1KM云图有2030×1354像素，分辨率更高的VIIRS则有超过6000×6000像素，像这种分辨率很高的极轨云图如果不经过裁剪，渲染时间会较长，而且出图后也很难上传到论坛/贴吧等平台。裁剪数据会用到python的“切片”功能，注意这里需要用到多维（二维）数组的切片，而不是普通的一维切片。
以bt31为例，二维数组的切片方法如下：
bt31[a:b, c:d]
a，b代表第一个方向上的起止序号；c，d代表第二个方向，这样能裁剪出以TC为中心的近似平行四边形的图像，我现在画的大部分云图也是这种。
如果需要得到长方形图像可以直接通过经纬度裁剪，未使用Cartopy的话可以直接用plt.xlim，plt.ylim输入经纬度进行裁剪，如果使用了Cartopy可以使用set_extent，若直接通过经纬度裁剪figsize也需要相应地调整。
2. 眼温读取
这个没什么技术难度，通过适当的切片使TC眼区为图像中亮温最高的部分，随后使用np.max读取即可。
3. 海岸线/不同投影
一般会选择Cartopy实现海岸线绘制和等经纬以外的投影。
import cartopy.crs as ccrs
import cartopy.feature as cfeature
fig = plt.figure(figsize=(londiff, latdiff), dpi=150) #给创建的画布一个名字（fig）
ax = fig.subplots(1, 1, subplot_kw={'projection': ccrs.Mercator()}) #设置子图为墨卡托投影，使高纬TC形态不被“压扁”
ax.add_feature(cfeature.COASTLINE.with_scale('10m'), lw=1) #添加海岸线
ax.pcolormesh(newlon, newlat, bt31, vmin=-100, vmax=50, cmap=LinearSegmentedColormap('1', bddata), transform=ccrs.PlateCarree()) #加上transform参数，明确数据原先是在等经纬网格上定义的

同一张Patricia，使用了墨卡托投影并加了墨西哥海岸线（不过画非高纬TC我平时还是最习惯等经纬）