前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >Python绘制气象实用地图(附代码和测试数据)

Python绘制气象实用地图(附代码和测试数据)

作者头像
MeteoAI
发布2019-08-12 14:41:33
14.8K3
发布2019-08-12 14:41:33
举报
文章被收录于专栏:MeteoAIMeteoAI

前面的推文对于常用的Python绘图工具都有了一些介绍,在这里就不赘述了。本文主要就以下几个方面:“中国区域绘图”“包含南海”“兰伯特投影带经纬度标签”“基于salem的mask方法”“进阶中国区域mask方法”“进阶省份mask方法”。对日常的实用需求能够在一定程度上满足。后续就Python在气象常用的统计方法(显著性检验)、合成分析、多变量叠加绘图再进行推送,敬请期待!

简单粗暴,Just show you my code!,细节暂不做过多分析,有问题可以探讨。数据、中文字体、地图shapefile文件、代码后文全部提供。使用建议,根据提示把缺失的库使用pip install xxx /conda install xxx /python setup.py install;安装完备,Python环境管理只推荐conda来统一管理。IDE推荐:PyCharm(有教育版)本地/服务器远程、Jupyter notebook。


•绘制兰勃脱投影的中国区域(包含南海子图)

•Mask掉海洋部分的兰勃脱投影(包含南海子图)

•基于salem的白化

•中国区域白化(包含南海子图)

•单独省份区域白化


1.绘制兰勃脱投影的中国区域(包含南海子图):

代码语言:javascript
复制
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from copy import copy
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom



def find_side(ls, side):
    """
 Given a shapely LineString which is assumed to be rectangular, return the
 line corresponding to a given side of the rectangle.

 """
    minx, miny, maxx, maxy = ls.bounds
    points = {'left': [(minx, miny), (minx, maxy)],
              'right': [(maxx, miny), (maxx, maxy)],
              'bottom': [(minx, miny), (maxx, miny)],
              'top': [(minx, maxy), (maxx, maxy)],}
    return sgeom.LineString(points[side])

def lambert_xticks(ax, ticks):
    """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    te = lambda xy: xy[0]
    lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
    xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
    ax.xaxis.tick_bottom()
    ax.set_xticks(xticks)
    ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])

def lambert_yticks(ax, ticks):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
    ax.yaxis.tick_left()
    ax.set_yticks(yticks)
    ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])

def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    """Get the tick locations and labels for an axis of a Lambert Conformal projection."""
    outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
    axis = find_side(outline_patch, tick_location)
    n_steps = 30
    extent = ax.get_extent(ccrs.PlateCarree())
    _ticks = []
    for t in ticks:
        xy = line_constructor(t, n_steps, extent)
        proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
        xyt = proj_xyz[..., :2]
        ls = sgeom.LineString(xyt.tolist())
        locs = axis.intersection(ls)
        if not locs:
            tick = [None]
        else:
            tick = tick_extractor(locs.xy)
        _ticks.append(tick[0])
    # Remove ticks that aren't visible: 
    ticklabels = copy(ticks)
    while True:
        try:
            index = _ticks.index(None)
        except ValueError:
            break
        _ticks.pop(index)
        ticklabels.pop(index)
    return _ticks, ticklabels


# Load the border data, CN-border-La.dat is downloaded from
# https://gmt-china.org/data/CN-border-La.dat
with open('CN-border-La.dat') as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]


# Set figure size
proj = ccrs.LambertConformal(central_longitude=105, central_latitude=90,
                             false_easting=400000, false_northing=400000)#,standard_parallels=(46, 49))

fig = plt.figure(figsize=[10, 8],frameon=True)


# Set projection and plot the main figure
ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=proj)
# Set figure extent
ax.set_extent([80, 130, 15, 55],crs=ccrs.PlateCarree())

# Plot border lines
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',
            transform=ccrs.Geodetic())

# Add ocean, land, rivers and lakes
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))


# *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 = np.arange(80,130,10)
# yticks = np.arange(15,55,5)
xticks = [55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155, 165]
yticks = [0 , 5 , 10, 15, 20, 25 , 30 , 35 , 40 , 45 , 50 , 55 , 60 , 65]
ax.gridlines(xlocs=xticks, ylocs=yticks,linestyle='--',lw=1,color='dimgrey')

# 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)


#Plot South China Sea as a subfigure
sub_ax = fig.add_axes([0.754, 0.107, 0.14, 0.155],
                      projection=ccrs.LambertConformal(central_latitude=90,
                                                       central_longitude=115))
# Add ocean, land, rivers and lakes
sub_ax.add_feature(cfeature.OCEAN.with_scale('50m'))
sub_ax.add_feature(cfeature.LAND.with_scale('50m'))
sub_ax.add_feature(cfeature.RIVERS.with_scale('50m'))
sub_ax.add_feature(cfeature.LAKES.with_scale('50m'))
# Plot border lines
for line in borders:
    sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
                transform=ccrs.Geodetic())
# Set figure extent
sub_ax.set_extent([105, 125, 0, 25],crs=ccrs.PlateCarree())

# Show figure
plt.savefig("China_nine_lines_lambert.png", dpi=500, bbox_inches='tight')
plt.show()

出图:

2.Mask掉海洋部分的兰勃脱投影(包含南海子图)

代码语言:javascript
复制
# https://gnss.help/2018/04/24/cartopy-gallery/index.html
# http://www.meteo.mcgill.ca/~cmccray/python.html
# https://www.jianshu.com/p/6506691f0788
# https://nbviewer.jupyter.org/gist/ajdawson/dd536f786741e987ae4e

import numpy as np
import os,sys
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from copy import copy
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom


def find_side(ls, side):
    """
 Given a shapely LineString which is assumed to be rectangular, return the
 line corresponding to a given side of the rectangle.

 """
    minx, miny, maxx, maxy = ls.bounds
    points = {'left': [(minx, miny), (minx, maxy)],
              'right': [(maxx, miny), (maxx, maxy)],
              'bottom': [(minx, miny), (maxx, miny)],
              'top': [(minx, maxy), (maxx, maxy)],}
    return sgeom.LineString(points[side])

def lambert_xticks(ax, ticks):
    """Draw ticks on the bottom x-axis of a Lambert Conformal projection."""
    te = lambda xy: xy[0]
    lc = lambda t, n, b: np.vstack((np.zeros(n) + t, np.linspace(b[2], b[3], n))).T
    xticks, xticklabels = _lambert_ticks(ax, ticks, 'bottom', lc, te)
    ax.xaxis.tick_bottom()
    ax.set_xticks(xticks)
    ax.set_xticklabels([ax.xaxis.get_major_formatter()(xtick) for xtick in xticklabels])

def lambert_yticks(ax, ticks):
    """Draw ricks on the left y-axis of a Lamber Conformal projection."""
    te = lambda xy: xy[1]
    lc = lambda t, n, b: np.vstack((np.linspace(b[0], b[1], n), np.zeros(n) + t)).T
    yticks, yticklabels = _lambert_ticks(ax, ticks, 'left', lc, te)
    ax.yaxis.tick_left()
    ax.set_yticks(yticks)
    ax.set_yticklabels([ax.yaxis.get_major_formatter()(ytick) for ytick in yticklabels])

def _lambert_ticks(ax, ticks, tick_location, line_constructor, tick_extractor):
    """Get the tick locations and labels for an axis of a Lambert Conformal projection."""
    outline_patch = sgeom.LineString(ax.outline_patch.get_path().vertices.tolist())
    axis = find_side(outline_patch, tick_location)
    n_steps = 30
    extent = ax.get_extent(ccrs.PlateCarree())
    _ticks = []
    for t in ticks:
        xy = line_constructor(t, n_steps, extent)
        proj_xyz = ax.projection.transform_points(ccrs.Geodetic(), xy[:, 0], xy[:, 1])
        xyt = proj_xyz[..., :2]
        ls = sgeom.LineString(xyt.tolist())
        locs = axis.intersection(ls)
        if not locs:
            tick = [None]
        else:
            tick = tick_extractor(locs.xy)
        _ticks.append(tick[0])
    # Remove ticks that aren't visible: 
    ticklabels = copy(ticks)
    while True:
        try:
            index = _ticks.index(None)
        except ValueError:
            break
        _ticks.pop(index)
        ticklabels.pop(index)
    return _ticks, ticklabels


# Load the border data, CN-border-La.dat is downloaded from
# https://gmt-china.org/data/CN-border-La.dat
with open('CN-border-La.dat') as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]


# Set figure size
proj = ccrs.LambertConformal(central_longitude=105, central_latitude=90,
                             false_easting=400000, false_northing=400000)#,standard_parallels=(46, 49))

fig = plt.figure(figsize=[10, 8],frameon=True)

# Set projection and plot the main figure
ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=proj)
# Set figure extent
ax.set_extent([80, 130, 15, 55],crs=ccrs.PlateCarree())
# ax.set_title('ChinaChinaChinaChinaChinaChinaChinaChinaChinaChina',fontsize=20)       
# Plot border lines
for line in borders:
    ax.plot(line[0::2], line[1::2], '-', lw=1.5, color='k',
            transform=ccrs.Geodetic())

# Add ocean, land, rivers and lakes
ax.add_feature(cfeature.OCEAN.with_scale('50m'))
ax.add_feature(cfeature.LAND.with_scale('50m'))
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))


# *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 = np.arange(80,130,10)
# yticks = np.arange(15,55,5)
xticks = [55, 65, 75, 85, 95, 105, 115, 125, 135, 145, 155, 165]
yticks = [0 , 5 , 10, 15, 20, 25 , 30 , 35 , 40 , 45 , 50 , 55 , 60 , 65]
ax.gridlines(xlocs=xticks, ylocs=yticks,linestyle='--',lw=1,color='dimgrey')

# 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)

# https://stackoverflow.com/questions/30030328/correct-placement-of-colorbar-relative-to-geo-axes-cartopy
# ax.set_aspect('auto', adjustable=None)

lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks)

#Plot South China Sea as a subfigure
sub_ax = fig.add_axes([0.592, 0.189, 0.14, 0.155],
                      projection=ccrs.LambertConformal(central_latitude=90,
                                                       central_longitude=115))
# Add ocean, land, rivers and lakes
sub_ax.add_feature(cfeature.OCEAN.with_scale('50m'))
sub_ax.add_feature(cfeature.LAND.with_scale('50m'))
sub_ax.add_feature(cfeature.RIVERS.with_scale('50m'))
sub_ax.add_feature(cfeature.LAKES.with_scale('50m'))
# Plot border lines
for line in borders:
    sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
                transform=ccrs.Geodetic())
# Set figure extent
sub_ax.set_extent([105, 125, 0, 25],crs=ccrs.PlateCarree())

def mask(ds, label='land'):
    landsea = xr.open_dataset('landsea.nc')
    landsea = landsea['LSMASK']
    # --ds和地形数据分辨率不一致,需将地形数据插值
    landsea = landsea.interp(lat=ds.latitude.values, lon=ds.longitude.values)
    # --利用地形掩盖海陆数据
    ds.coords['mask'] = (('latitude', 'longitude'), landsea.values)
    print(ds.mask)
    if label == 'land':
        ds = ds.where(ds.mask < 0.8) #可以尝试调整default:0.8
    elif label == 'ocean':
        ds = ds.where(ds.mask > 0.3) #可以尝试调整default:0.2
    return ds

ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
time = ds.time
temp = (ds['t2m'] - 273.15) # 把温度转换为℃
# 区域选择
lon_range = lon[(lon>60) & (lon<150)]
lat_range = lat[(lat>10) & (lat<60)]
temp_region = temp.sel(longitude=lon_range, latitude=lat_range, time='2018-02-01')
temp_mask = mask(temp_region, 'ocean')
# 设置colorbar 
#get size and extent of axes:

cbar_kwargs = {
    'orientation': 'vertical',
    'label': '2m temperature (℃)',
    'shrink': 0.8,
    'ticks': np.arange(-30, 30 + 5, 5),
    'pad': 0.05,
    'shrink': 0.65
}


# 画图
levels = np.arange(-30, 30 + 1, 1)
temp_mask.plot.contourf(ax=ax, levels=levels, cmap='Spectral_r',
                        cbar_kwargs=cbar_kwargs, transform=ccrs.PlateCarree())
# 设置标题的在代码中放置的位置很关键,注意不要放置在小图上或者新建画框了。
ax.set_title('China Map 2m Temperature',color='blue',fontsize= 20)

# Save & Show figure
# (filename, extension) = os.path.splitext(os.path.basename(sys.argv[0]))
plt.savefig("land_sea_mask_simple_lambert.png", dpi=500, bbox_inches='tight')
plt.show()

出图

3.基于salm的白化

代码语言:javascript
复制
import salem
import matplotlib.pyplot as plt
from salem.utils import get_demo_file
ds = salem.open_xr_dataset(get_demo_file('wrfout_d01.nc'))
t2 = ds.T2.isel(Time=2)
# t2.salem.quick_map()
t2_sub = t2.salem.subset(corners=((77., 20.), (97., 35.)), crs=salem.wgs84)
# t2_sub.salem.quick_map()
shdf = salem.read_shapefile(get_demo_file('world_borders.shp'))
shdf = shdf.loc[shdf['CNTRY_NAME'].isin(['Nepal', 'Bhutan'])]  # GeoPandas' GeoDataFrame
t2_sub = t2_sub.salem.subset(shape=shdf, margin=2)  # add 2 grid points
# t2_sub.salem.quick_map()
t2_roi = t2_sub.salem.roi(shape=shdf)
# t2_roi.salem.quick_map()

# Change the country borders
smap = t2_roi.salem.get_map(data=t2_roi-273.15, cmap='RdYlBu_r', vmin=-14, vmax=18)
_ = smap.set_topography(get_demo_file('himalaya.tif'))
smap.set_shapefile(shape=shdf, color='grey', linewidth=3)
smap.set_shapefile(countries=True, color='C3', linewidths=2)

# Add oceans and lakes
smap.set_shapefile(oceans=True)
smap.set_shapefile(rivers=True,linewidths=0.8,color='green')
smap.set_shapefile(lakes=True, facecolor='blue', edgecolor='blue')

# Change the lon-lat countour setting
smap.set_lonlat_contours(add_ytick_labels=True, interval=2, linewidths=0.5,
                         linestyles='--', colors='grey')

# Add a scalebar (experimental)
smap.set_scale_bar(location=(0.10, 0.10), add_bbox=False)

smap.set_points(91.1, 29.6)
smap.set_text(90.0, 30.1, 'Lhasa', fontsize=17)

smap.visualize()
plt.savefig("Himalayas_mask.png", dpi=500, bbox_inches='tight')
plt.show()

出图:

4.中国区域白化(包含南海子图)

代码语言:javascript
复制
#https://regionmask.readthedocs.io/en/stable/defined_landmask.html
#coding=utf-8
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
import cartopy.crs as ccrs
import shapefile
import xarray as xr
import matplotlib as mpl
from matplotlib.font_manager import FontProperties
import netCDF4 as nc
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import matplotlib
ZHfont = matplotlib.font_manager.FontProperties(fname='/Users/zhpfu/Documents/fonts/SimSun.ttf')

plt.rcParams.update({'font.size':20})


fig = plt.figure(figsize=[12,18]) 
ax = fig.add_subplot(111)



sf = shapefile.Reader("country1")
for shape_rec in sf.shapeRecords():
    if shape_rec.record[2] == 'China':#Hunan Sheng
        vertices = []
        codes = []
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)



def makedegreelabel(degreelist):
    labels=[str(x)+u'°E' for x in degreelist]
    return labels


ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃   [0,::-1,:]表示第一个时次、纬度反向

# ncdata=nc.Dataset('EC-Interim_monthly_2018.nc')
# data=ncdata.variables['t2m'][0,::-1,:] - 273.15 # 把温度转换为℃ 
# lat=ncdata.variables['latitude'][:]
# lon=ncdata.variables['longitude'][:]

nx=data.shape[1]
ny=data.shape[0]

m = Basemap(llcrnrlon=72.0,
    llcrnrlat=10.0,
    urcrnrlon=137.0,
    urcrnrlat=55.0,
    resolution = None, 
    projection = 'cyl')


cbar_kwargs = {
    'orientation': 'horizontal',
    'label': 'Temperature (℃)',
    'shrink': 0.02,
    'ticks': np.arange(-25, 25 + 1, 5),
    'pad': -0.28,
    'shrink': 0.95
}


# 画图
levels = np.arange(-25, 25 + 1, 1)   
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
# 设置标题的在代码中放置的位置很关键,注意不要放置在小图上或者新建画框了。

m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2)

for contour in cs.collections:
        contour.set_clip_path(clip)


# draw parallels and meridians.
# label parallels on right and top
# meridians on bottom and left
# http://code.activestate.com/recipes/578399-an-alternative-way-to-draw-parallels-and-meridians/
parallels = np.arange(10,55,10)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])  # ha= 'right'
meridians = np.arange(70,135,10)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])


#http://xarray.pydata.org/en/stable/plotting.html
plt.ylabel('')    #Remove the defult  lat / lon label  
plt.xlabel('')        


plt.rcParams.update({'font.size':30})
ax.set_title(u' 中国区域2018年1月平均气温',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature


plt.rcParams.update({'font.size':20})

with open('CN-border-La.dat') as src:
    context = src.read()
    blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
    borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]

#Plot South China Sea as a subfigure
sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],
                      projection=ccrs.LambertConformal(central_latitude=90,
                                                       central_longitude=115))
# Plot border lines
for line in borders:
    sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',
                transform=ccrs.Geodetic())
# Set figure extent
sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree())        


# Save & Show figure
plt.savefig("China_mask_T2m.png", dpi=300, bbox_inches='tight')
plt.show()

出图:

5.单独省份区域白化

代码语言:javascript
复制
#https://basemaptutorial.readthedocs.io/en/latest/clip.html
#coding=utf-8
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.pyplot as plt
from osgeo import gdal
import numpy as np
import cartopy.crs as ccrs
import shapefile
import matplotlib as mpl
import xarray as xr
from matplotlib.font_manager import FontProperties
import netCDF4 as nc
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import matplotlib
ZHfont = matplotlib.font_manager.FontProperties(fname='/Users/zhpfu/Documents/fonts/SimSun.ttf')


plt.rcParams.update({'font.size': 20})

# plt.clf()
fig = plt.figure(figsize=[12,18]) 

# fig = plt.figure(figsize=(16,8.5),dpi=300)# figsize=(16,9.6)   
ax = fig.add_subplot(111)


# set the projection to PlateCarree
# proj = ccrs.PlateCarree()
# ax = fig.add_subplot(1, 1, 1,projection = proj)

# ax=fig.add_axes([.2,.2,.6,.68]) #[*left*, *bottom*, *width*,*height*]

sf = shapefile.Reader("national_provinces_shp/hunan")
for shape_rec in sf.shapeRecords():
    if shape_rec.record[2] == 'Hunan Sheng':#Hunan Sheng
        vertices = []
        codes = []
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                vertices.append((pts[j][0], pts[j][1]))
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform=ax.transData)



def makedegreelabel(degreelist):
    labels=[str(x)+u'°E' for x in degreelist]
    return labels


ds = xr.open_dataset('EC-Interim_monthly_2018.nc')
lat = ds.latitude
lon = ds.longitude
data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃

nx=data.shape[1]
ny=data.shape[0]


# Hunan Province region = [108.6,114.3,24.5,30.2] # [west,east,south,north]

m = Basemap(llcrnrlon=108.6,
    llcrnrlat=24.5,
    urcrnrlon=114.3,
    urcrnrlat=30.2,
    resolution = None, 
    projection = 'cyl')


cbar_kwargs = {
    'orientation': 'horizontal',
    'label': 'Temperature (℃)',
    'shrink': 0.02,
    'ticks': np.arange(0, 10 + 1, 1),
    'pad': -0.01,
    'shrink': 0.95
}



# 画图
levels = np.arange(0, 8 + 1, 0.2)   
cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')
# 设置标题的在代码中放置的位置很关键,注意不要放置在小图上或者新建画框了。

m.readshapefile('national_provinces_shp/hunan','Hunan Map',color='k',linewidth=1.2)

for contour in cs.collections:
        contour.set_clip_path(clip)

# draw parallels and meridians.
# label parallels on right and top
# meridians on bottom and left
# http://code.activestate.com/recipes/578399-an-alternative-way-to-draw-parallels-and-meridians/
parallels = np.arange(24.5,30.2,1.0)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3])  # ha= 'right'
meridians = np.arange(108.6,114.3,1.0)
m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3])


#http://xarray.pydata.org/en/stable/plotting.html
plt.ylabel('')    #Remove the defult  lat / lon label  
plt.xlabel('')


plt.rcParams.update({'font.size':30})
ax.set_title(u' 湖南省2018年1月平均气温',color='blue',fontsize= 25 ,fontproperties=ZHfont) # 2m Temperature

bill0 = 111.69
tip0  =  29.05
plt.scatter(bill0, tip0,marker='.',s=100 ,color ="blue")


bill2 =  113.0
tip2  = 28.21
plt.scatter(bill2, tip2,marker='*',s=150 ,color ="orange" )

plt.rcParams.update({'font.size':30})
plt.text(bill0-0.4, tip0+0.2, u"常德市", fontsize= 20 ,fontproperties=ZHfont)
plt.text(bill2-0.4, tip2+0.2, u"长沙市", fontsize= 20 ,fontproperties=ZHfont)

# Save & Show figure
plt.savefig("Hunan_mask_T2m.png", dpi=300, bbox_inches='tight')
plt.show()

出图:

(老家在湖南,?,优先示例!)

测试数据和代码:链接:https://pan.baidu.com/s/18R6RWYhi5p_wMbMrdKzw2g 密码:jwil

References

  • https://basemaptutorial.readthedocs.io/en/latest/clip.html
  • http://code.activestate.com/recipes/578399-an-alternative-way-to-draw-parallels-and-meridians/
  • https://stackoverflow.com/questions/30030328/correct-placement-of-colorbar-relative-to-geo-axes-cartopy
  • https://gnss.help/2018/04/24/cartopy-gallery/index.html
  • http://www.meteo.mcgill.ca/~cmccray/python.html
  • https://www.jianshu.com/p/6506691f0788
  • https://nbviewer.jupyter.org/gist/ajdawson/dd536f786741e987ae4e
  • https://gmt-china.org/data/CN-border-La.dat

本文同步更新在微信公众号『气象学家』

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2019-08-10,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 MeteoAI 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1.绘制兰勃脱投影的中国区域(包含南海子图):
  • 2.Mask掉海洋部分的兰勃脱投影(包含南海子图)
  • 3.基于salm的白化
  • 4.中国区域白化(包含南海子图)
  • 5.单独省份区域白化
  • References
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档