专栏首页气象杂货铺绘图|解决Cartopy Lambert投影坐标轴标签设置问题

绘图|解决Cartopy Lambert投影坐标轴标签设置问题

python中有两个使用最频繁的地图绘图库:Basemap和Cartopy,两者各有优劣。由于Cartopy和matplotlib的兼容性更好,并且用户友好度更高,开始逐渐被人接受。但是Cartopy也有一些缺点,其中之一就是在设置坐标轴标签的时候对于非矩形投影无法设置标签,比如Lambert投影。

对于不受投影限制的绘图可以转换为PlateCarree投影或者Mercator投影,但对于有投影限制的绘图,比如WRF模式的后处理(虽然WRF模式也支持Mercator投影,但是大多数情况下还是使用的Lambert投影)就受限了。

在互联网游荡的时候偶然发现了一个用于解决此问题的脚本[注1],然后测试了一下,发现基本能够完美解决Cartopy Lambert投影标签设置的问题。以下是脚本:

from copy import copy
import numpy as np
import shapely.geometry as sgeom
import cartopy.crs as ccrs

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 ticks 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

以下是测试结果:

import cartopy.crs as ccrs
import cartopy.feature as cfeature  
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import matplotlib.pyplot as plt

proj = ccrs.LambertConformal(central_longitude=117, central_latitude=26,
                             standard_parallels=(30, 60))

shpfn = 'bou2_4l.shp'
reader = shpreader.Reader(shpfn)
states_provinces = cfeature.ShapelyFeature(reader.geometries(), crs = ccrs.PlateCarree(), edgecolor = 'face', facecolor = 'None')

fig = plt.figure(figsize=(9, 9), frameon=True)

ax = fig.add_axes([0.08, 0.05, 0.8, 0.94], projection=proj)

# Label axes of a Mercator projection without degree symbols in the labels
# and formatting labels to include 1 decimal place:
ax.set_extent([85, 135, 0, 55], crs=ccrs.PlateCarree())

#    ax.coastlines(resolution='50m')
ax.add_feature(states_provinces, edgecolor='black')
# *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 = list(range(60, 176, 10))
yticks = list(range(-5, 66, 10))
ax.gridlines(xlocs=xticks, ylocs=yticks)

# 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)
# xticks and yticks only is list
lambert_xticks(ax, xticks)  
lambert_yticks(ax, yticks)

plt.show()

在设置坐标轴标签时仍然会存在一些小问题,但是这些都可以通过更改设置解决。上述提到的方法能够解决标签标注的问题,但是对numpy的支持不是很好,但是只需要进行一定的更改即可。


注1: https://gist.github.com/ajdawson/dd536f786741e987ae4e

注2: https://github.com/bugsuse/DataAnalysis/blob/master/Cartopy%E4%B8%ADLambert%E6%8A%95%E5%BD%B1%E8%AE%BE%E7%BD%AEticks.ipynb

本文分享自微信公众号 - 气象杂货铺(meteogs),作者:bugsuse

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-10-05

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 建议收藏!Matplotlib常见组件设置整理

    继上一篇文章为大家介绍了plt和ax绘图的区别后,这篇文章结合我自己的一些使用经历,为大家整理了Matplotlib中比较常用的一些组件设置。

    bugsuse
  • 使用 Cartopy 和 netCDF4 可视化 WRF 模式数据

    对比使用 Basemap,gdal 和 Cartopy,netCDF4 读取 WRF 模式数据并绘图。

    bugsuse
  • 真・WRF模式后处理之Python版

    WRF模式是数值天气预报和大气模拟系统,其开发目的就是用语研究和实际应用。运行WRF模式时,可以利用多种初始场数据来驱动,然后配置好选项之后便可以模拟天气过程(...

    bugsuse
  • 全文检索Solr集成HanLP中文分词

    以前发布过HanLP的Lucene插件,后来很多人跟我说其实Solr更流行(反正我是觉得既然Solr是Lucene的子项目,那么稍微改改配置就能支持Solr),...

    IT小白龙
  • JBPM4.4(1)-简单工程的搭建

    源码下载 https://anonsvn.jboss.org/repos/jbpm/jbpm4/ JBPM是什么? jBPM是一个可扩展、灵活的流程引擎, 它可...

    cloudskyme
  • python3装饰器

    装饰器本质其实就是一个函数, 可以让其它函数不改动源代码的情况下增加其他新功能, 比如网站经常需要的权限校验等场景

    py3study
  • Linux命令之time——计算命令运行时间

    linux下time命令可以获取到一个程序的执行时间,包括程序的实际运行时间(real time),以及程序运行在用户态的时间(user time)和内核态的时...

    浩Coding
  • JS实现网页简体繁体字转换功能

    在网页中经常会遇到将简体字转换成繁体字,方便于其他同胞查看。网页中实现简体中文转换成繁体字方法,今天分享给大家,此方法借鉴于他人博客;

    申霖
  • java中集合关系图

    黑泽君
  • FPGA设计心得(9)基于DDS IP核的任意波形发生器设计

    数据手册[1]博客首页[2]花了几个小时写了这篇博客,不得不说的是了解的还是皮毛而已,但尽力写的详细点,这比较适合新手,老手可以忽略繁琐的部分。注:学习交流使用...

    Reborn Lee

扫码关注云+社区

领取腾讯云代金券