Python气象处理绘图第一弹--白化海洋部分

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

---

前言

我是学习python处理气象数据以及绘图的小白。以前常用Matlab处理数据，利用ncl进行绘图。但是MATLAB对于批量nc文件的处理很慢，比较消耗时间，于是我转战python啦。

一、白化是什么？

在气象绘图过程中，我们常常需要对不感兴趣的部分进行遮掩。于是需要对不感兴趣的部分进行白化处理。目前比较常用的是maskout这一由平流层的萝卜首发、晋陵小生优化的 maskout 模块白化功能。这一功能对于区域的 白化非常有用，但我此次仅对海洋数据进行白化，我就参考了公众号好奇心 Log 推送的 salem 库包下的白化。

二、使用步骤

首先需要对salem程序包进行安装，详情可以度娘解决。

1.引入库

代码如下（示例）：

# -*- coding: utf-8 -*-
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.mpl.ticker as cticker
from scipy import stats
from matplotlib.font_manager import FontProperties
import warnings
#利用salem 进行陆地掩膜
import salem
import fiona
from cartopy.io.shapereader import Reader, natural_earth
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import geopandas

2.读入数据

代码如下（示例）：

shp_path = 'F:/mapshp/landmask/'
shp = geopandas.read_file(shp_path + 'ne_10m_land_scale_rank.shp')
mf1=xr.open_mfdataset('control/*.nc')
#修改nc文件中经度0-360~-180-180，要使用salem进行白化必须使用标准的经纬度
long_name = 'lon'  #你的nc文件中经度的命名
mf1['longitude_adjusted'] = xr.where(
    mf1[long_name] > 180,
    mf1[long_name] - 360,
    mf1[long_name])
mf1 = (
    mf1
    .swap_dims({long_name: 'longitude_adjusted'})
    .sel(**{'longitude_adjusted': sorted(mf1.longitude_adjusted)})
    .drop(long_name))
mf1 = mf1.rename({'longitude_adjusted': long_name})
#提取经纬度
lat = mf1['lat']
lon = mf1['lon']
tsc=mf1['TS'] - 273.15
tsc_clim=tsc.mean(dim=['time'])
tsc_climmask = tsc_clim.salem.roi(shape=shp)
def plot_with_map2(axes,img_extent,spec1,spec2,title):
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    axes.set_extent(img_extent, crs=ccrs.PlateCarree())
    axes.add_feature(cfeature.COASTLINE.with_scale('50m')) 
    axes.set_xticks(np.arange(leftlon,rightlon+spec1,spec1), crs=ccrs.PlateCarree())
    axes.set_yticks(np.arange(lowerlat,upperlat+spec2,spec2), crs=ccrs.PlateCarree())
    lon_formatter = cticker.LongitudeFormatter()
    lat_formatter = cticker.LatitudeFormatter()
    axes.xaxis.set_major_formatter(lon_formatter)
    axes.yaxis.set_major_formatter(lat_formatter)
    axes.set_title(title,fontsize=14,loc='left')
fig1 = plt.figure(figsize=(16,8))
proj = ccrs.PlateCarree(central_longitude=0)
leftlon, rightlon, lowerlat, upperlat = (-180,180,-90,90)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
ax1=plt.subplot(1,1,1,projection=proj)
tsfig = ax1.contourf(lon,lat, tsc_climmask,levels =np.arange(-60,30,0.5), 
                      extend = 'both', transform=ccrs.PlateCarree(), cmap=plt.cm.RdBu_r)
plot_with_map2(ax1,img_extent,60,30,'TS_C')

---

总结

在使用salem库的过程是使不感兴趣的格点值等于nan
用salem.roi提取区域的时候，维度名称要用标准的lat，lon，如果出现lat_0，lon_0是无法识别的。
注意对与geopandas包的安装