{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# NCL_overlay_11b.py\nThis script illustrates the following concepts:\n - Overlaying vectors and filled contours on a map\n - Masking out particular areas in a map\n - Subsetting a color map\n\nSee following URLs to see the reproduced NCL plot & script:\n - Original NCL script: https://www.ncl.ucar.edu/Applications/Scripts/overlay_11.ncl\n - Original NCL plot: https://www.ncl.ucar.edu/Applications/Images/overlay_11_lg.png\n\nThis script shows how to overlay contours and vectors on a map,\nbut with the contours limited to specific areas, and the vectors\nnot limited.\n\nThe point of this script is to show how to mask contours against\na geographical boundary, but in a way that allows them to be drawn\nup to the boundary location. This is unlike the shapefile masking\nexamples, where grid points are set to missing if they fall\noutside a boundary, and hence you can get blocky features close\nto the boundary.\n\nWith Python's matplotlib, there are 2 general approaches to\naccomplishing this:\n\na. You can \"cover\" (i.e., \"over lay\") geographical features on\n top of other plots using the ``zorder`` parameter to most\n rendered objects.\n\nb. You can \"clip\" a plot object with a geographical boundary.\n\nThis example demonstrates approach (b).\n\nSee NCL_overlay_11a.py for demonstration of approach (a).\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import packages:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import xarray as xr\nimport numpy as np\n\nfrom matplotlib import pyplot as plt\nfrom matplotlib.patches import PathPatch\n\nfrom cartopy.feature import ShapelyFeature, OCEAN, LAKES, LAND\nfrom cartopy.crs import PlateCarree\nfrom cartopy.mpl.patch import geos_to_path\nfrom cartopy.io.shapereader import Reader as ShapeReader, natural_earth\n\nimport geocat.datafiles as gdf\nfrom geocat.viz import cmaps as gvcmaps\nfrom geocat.viz import util as gvutil" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read in data:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Open a netCDF data file using xarray default engine and load the data into xarrays, as well as extract slices for\n# ``time=0`` and the ``lev=500`` hPa level\nds = xr.open_dataset(gdf.get(\"netcdf_files/uvt.nc\")).sel(time=0, lev=500)\n\n# For convenience only, extract the U,V,T and lat and lon variables\nU = ds[\"U\"]\nV = ds[\"V\"]\nT = ds[\"T\"]\n\nlat = ds[\"lat\"]\nlon = ds[\"lon\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Construct shape boundaries\n\nUsing Cartopy's interface to the Natural Earth Collection of shapefiles\nand geographical shape data, we construct the geographical boundaries\nthat we are interested in displaying, namely the country borders of China\nand Taiwan, the borders of Chinese provinces, and all land borders *without*\nChina or Taiwan.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Download the Natural Earth shapefile for country boundaries at 10m resolution\nshapefile = natural_earth(category='cultural',\n resolution='10m',\n name='admin_0_countries')\n\n# Sort the geometries in the shapefile into Chinese/Taiwanese or other\ncountry_geos = []\nother_land_geos = []\nfor record in ShapeReader(shapefile).records():\n if record.attributes['ADMIN'] in ['China', 'Taiwan']:\n country_geos.append(record.geometry)\n else:\n other_land_geos.append(record.geometry)\n\n# Define map projection to allow Cartopy to transform ``lat`` and ``lon`` values accurately into points on the\n# matplotlib plot canvas.\nprojection = PlateCarree()\n\n# Define a Cartopy Feature for the country borders and the land mask (i.e.,\n# all other land) from the shapefile geometries, so they can be easily plotted\ncountries = ShapelyFeature(country_geos,\n crs=projection,\n facecolor='none',\n edgecolor='black',\n lw=1.5)\nland_mask = ShapelyFeature(other_land_geos,\n crs=projection,\n facecolor='white',\n edgecolor='none')\n\n# Download the Natural Earth shapefile for the states/provinces at 10m resolution\nshapefile = natural_earth(category='cultural',\n resolution='10m',\n name='admin_1_states_provinces')\n\n# Extract the Chinese province borders\nprovince_geos = [\n record.geometry\n for record in ShapeReader(shapefile).records()\n if record.attributes['admin'] == 'China'\n]\n\n# Define a Cartopy Feature for the province borders, so they can be easily plotted\nprovinces = ShapelyFeature(province_geos,\n crs=projection,\n facecolor='none',\n edgecolor='black',\n lw=0.25)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Generate figure (set its size (width, height) in inches) and axes using Cartopy\nfig = plt.figure(figsize=(10, 10))\nax = plt.axes(projection=projection)\n\nax.set_extent([100, 145, 15, 55], crs=projection)\n\n# Draw the ocean and lake features\nax.add_feature(OCEAN.with_scale('50m'), edgecolor='black', lw=1)\nax.add_feature(LAKES.with_scale('50m'), edgecolor='black', lw=1)\n\n# Define the contour levels (T)\nclevs = np.arange(228, 273, 4, dtype=float)\n\n# Import an NCL colormap, truncating it by using geocat.viz.util convenience function\nnewcmp = gvutil.truncate_colormap(gvcmaps.BkBlAqGrYeOrReViWh200,\n minval=0.1,\n maxval=0.6,\n n=len(clevs))\n\n# Draw the contour plot, \"clipped\" to the country boundaries\n# (NOTE: There are multiple closed polygons representing the boundaries of the\n# countries. This is both because there are 2 country borders being used\n# to clip the contour plot, but also because China consists of many islands.\n# As a result, we have to loop over *all closed paths* and construct a\n# matplotlib patch object that we can use the clip the contour plot.)\nfor path in geos_to_path(country_geos):\n patch = PathPatch(path,\n transform=ax.transData,\n facecolor='none',\n edgecolor='black',\n lw=1.5)\n\n # Draw the patch on the plot\n ax.add_patch(patch)\n\n # Draw the contour plot\n # (NOTE: Because this line is in the loop over closed paths, the contour plot\n # is being drawn for each closed path. This has to be done because\n # matplotlib cannot apply *multiple* closed paths at the same time to\n # to the same plot. Hence, for each closed path, we need to generate\n # another contour plot and clip that contour plot with the patch. In\n # other words, every island on this plot corresponds to its own\n # contour plot!)\n cf = ax.contourf(lon, lat, T, levels=clevs, cmap=newcmp)\n\n # Clip each contour of the contour plot\n # (NOTE: Each contour of the contour plot is actually its own \"plot\". There\n # is no easy mechanism in matplotlib to clip the entire contour plot\n # at once, so we must loop through the \"collections\" in the contour\n # plot and clip each one separately.)\n for col in cf.collections:\n col.set_clip_path(patch)\n\n# Add horizontal colorbar\ncax = plt.axes((0.14, 0.08, 0.74, 0.02))\ncbar = plt.colorbar(cf,\n ax=ax,\n cax=cax,\n ticks=clevs[1:-1],\n drawedges=True,\n orientation='horizontal')\ncbar.ax.tick_params(labelsize=12)\n\n# Draw the province borders\nax.add_feature(provinces)\n\n# Draw the quiver plot (and its key)\nQ = ax.quiver(lon,\n lat,\n U,\n V,\n color='black',\n width=.003,\n scale=600.,\n headwidth=3.75)\nrect = plt.Rectangle((142, 52),\n 3,\n 3,\n facecolor='mediumorchid',\n edgecolor=None,\n zorder=1)\nax.add_patch(rect)\nax.quiverkey(Q,\n 0.9675,\n 0.95,\n 30,\n '30',\n labelpos='N',\n color='black',\n coordinates='axes',\n fontproperties={'size': 14},\n labelsep=0.1)\n\n# Draw the '500hPa' label at the top left of the plot\nprops = dict(facecolor='white', edgecolor='none', alpha=0.8)\nax.text(105,\n 52.7,\n '500hPa',\n transform=projection,\n fontsize=18,\n ha='center',\n va='center',\n color='mediumorchid',\n bbox=props)\n\n# Use geocat.viz.util convenience function to set axes tick values\ngvutil.set_axes_limits_and_ticks(ax,\n xticks=[100, 120, 140],\n yticks=[20, 30, 40, 50])\n\n# Use geocat.viz.util convenience function to make plots look like NCL plots by using latitude, longitude tick labels\ngvutil.add_lat_lon_ticklabels(ax)\n\n# Use geocat.viz.util convenience function to add minor and major tick lines\ngvutil.add_major_minor_ticks(ax,\n x_minor_per_major=4,\n y_minor_per_major=5,\n labelsize=18)\n\n# Use geocat.viz.util convenience function to add main title as well as titles to left and right of the plot axes.\ngvutil.set_titles_and_labels(ax,\n lefttitle=\"Temp\",\n lefttitlefontsize=20,\n righttitle=\"Wind\",\n righttitlefontsize=20)\n\n# Show the plot\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }