01_example-CCMP-notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Access CCMP data on Pangeo\n",
    "\n",
    "This tutorial shows how to access the CCMP wind speed data, makes a few plots of the data, and then compares it to a buoy wind speed.\n",
    "\n",
    "## CCMP Winds in a cloud-optimized-format for Pangeo\n",
    "\n",
    "The Cross-Calibrated Multi-Platform (CCMP) Ocean Surface Wind Vector Analyses is part of the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) Program. MEaSUREs, develops consistent global- and continental-scale Earth System Data Records by supporting projects that produce data using proven algorithms and input.  If you use this data, please give [credit](https://podaac.jpl.nasa.gov/MEaSUREs-CCMP?sections=about).  For more information, please review the [documentation](https://podaac-tools.jpl.nasa.gov/drive/files/allData/ccmp/L2.5/docs/ccmp_users_guide.pdf). Please note that this data is not recommended for trend calculations.\n",
    "\n",
    "Data is NRT from 4/1/2019 - 10/4/2020.\n",
    "The 6-hourly VAM analyses are centered at 0,6,12 and 18z\n",
    "\n",
    "## Accessing cloud satellite data\n",
    "\n",
    "- CCMP zarr conversion funding: Interagency Implementation and Advanced Concepts Team [IMPACT](https://earthdata.nasa.gov/esds/impact) for the Earth Science Data Systems (ESDS) program and AWS Public Dataset Program\n",
    "  \n",
    "### Credit:\n",
    "- [Chelle Gentemann](mailto:cgentemann@faralloninstitute.org), [Farallon Institute](http://www.faralloninstitute.org/), [Twitter](https://twitter.com/ChelleGentemann) - creation of Zarr data store and tutorial\n",
    "- [Charles Blackmon Luca](mailto:blackmon@ldeo.columbia.edu), [LDEO](https://www.ldeo.columbia.edu/) - help with moving to Pangeo storage and intake update\n",
    "- [Willi Rath](mailto:wrath@geomar.de), [GEOMAR](https://www.geomar.de/en/), [Twitter](https://twitter.com/RathWilli) - motivated CG to move data to Pangeo!\n",
    "- [Tom Augspurger](mailto:taugspurger@microsoft.com), [Microsoft](https://www.linkedin.com/in/tom-augspurger-01269315/),[Twitter](https://twitter.com/TomAugspurger) - help updating the data.\n",
    "\n",
    "#### Zarr data format\n",
    "\n",
    " [Zarr](https://zarr.readthedocs.io/en/stable/)\n",
    "\n",
    "#### Data proximate computing\n",
    "These are BIG datasets that you can analyze on the cloud without downloading the data. You can run this on your phone, a Raspberry Pi, laptop, or desktop.   \n",
    "By using public cloud data, your science is reproducible and easily shared!\n",
    "\n",
    "#### To run this notebook\n",
    "\n",
    "Code is in the cells that have <span style=\"color: blue;\">In [  ]:</span> to the left of the cell and have a colored background\n",
    "\n",
    "To run the code:\n",
    "- option 1) click anywhere in the cell, then hold `shift` down and press `Enter`\n",
    "- option 2) click on the Run button at the top of the page in the dashboard\n",
    "\n",
    "Remember:\n",
    "- to insert a new cell below press `Esc` then `b`\n",
    "- to delete a cell press `Esc` then `dd`\n",
    "\n",
    "#### First start by importing libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#libs for reading data\n",
    "import xarray as xr\n",
    "import intake\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.dates import DateFormatter\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "#lib for dask gateway\n",
    "from dask_gateway import Gateway\n",
    "from dask.distributed import Client\n",
    "from dask import delayed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Start a cluster, a group of computers that will work together.\n",
    "\n",
    "(A cluster is the key to big data analysis on on Cloud.)\n",
    "\n",
    "- This will set up a [dask kubernetes](https://docs.dask.org/en/latest/setup/kubernetes.html) cluster for your analysis and give you a path that you can paste into the top of the Dask dashboard to visualize parts of your cluster.  \n",
    "- You don't need to paste the link below into the Dask dashboard for this to work, but it will help you visualize progress.\n",
    "- Try 20 workers to start (during the tutorial) but you can increase to speed things up later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gateway = Gateway()\n",
    "cluster = gateway.new_cluster()\n",
    "cluster.adapt(minimum=1, maximum=75)\n",
    "client = Client(cluster)\n",
    "cluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** ☝️ Don’t forget to click the link above or copy it to the Dask dashboard (orange icon) on the left to view the scheduler dashboard! **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialize Dataset\n",
    "\n",
    "Here we load the dataset from the zarr store. Note that this very large dataset (273 GB) initializes nearly instantly, and we can see the full list of variables and coordinates.\n",
    "\n",
    "### Examine Metadata\n",
    "\n",
    "For those unfamiliar with this dataset, the variable metadata is very helpful for understanding what the variables actually represent\n",
    "Printing the dataset will show you the dimensions, coordinates, and data variables with clickable icons at the end that show more metadata and size."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "cat_pangeo = intake.open_catalog(\"https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/master.yaml\")\n",
    "\n",
    "ds_ccmp = cat_pangeo.atmosphere.nasa_ccmp_wind_vectors.to_dask()\n",
    "\n",
    "#calculate wind speed and add attributes to new variable\n",
    "ds_ccmp['wspd'] = np.sqrt(ds_ccmp.uwnd**2 + ds_ccmp.vwnd**2)\n",
    "ds_ccmp.wspd.attrs=ds_ccmp.uwnd.attrs\n",
    "ds_ccmp.wspd.attrs['long_name']='wind speed at 10 meters'\n",
    "\n",
    "ds_ccmp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot a global image of the data on 7/4/2020 (randomly chosen day)\n",
    "\n",
    "``xarray`` makes plotting the data very easy.  A nice overview of plotting with xarray is [here](http://xarray.pydata.org/en/stable/plotting.html).  Details on [.plot](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.plot.html#xarray.DataArray.plot)\n",
    "\n",
    "Below we are plotting the variable `wspd` (wind speed) and `nobs` (number of observations) to show how there is data over the ocean but none over land. CCMP data includes 25 km ERA5 data over land."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "day = ds_ccmp.sel(time='2020-09-08T18')\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(18,6))\n",
    "\n",
    "day.wspd.plot(ax=ax[0])\n",
    "\n",
    "day.nobs.plot(ax=ax[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "ortho = ccrs.Orthographic(-110, 20)           # define target coordinate frame\n",
    "geo = ccrs.PlateCarree()                     # define origin coordinate frame\n",
    "\n",
    "plt.figure(figsize=(7,5))                    #set the figure size\n",
    "ax = plt.subplot(1, 1, 1, projection=ortho)  #create the axis for plotting\n",
    "q = day.uwnd.plot(ax=ax, \n",
    "                            transform = geo, \n",
    "                            cmap='RdBu_r', #'cmo.thermal', \n",
    "                            vmin=-5, \n",
    "                            vmax=5,cbar_kwargs={'label':'Eastward Wind Speed (m/s)'}) # plot a colormap in transformed coordinates cmap='OrRd', \n",
    "global_extent = ax.get_extent(crs=ccrs.PlateCarree())\n",
    "#ax.add_feature(cartopy.feature.LAND,color='gray')\n",
    "ax.add_feature(cartopy.feature.COASTLINE)\n",
    "\n",
    "#plt.savefig('./../../figures/ccmp_wind_example.png')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot a timeseries of the average wind speed over a region\n",
    "\n",
    "- use [.sel](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.sel.html#xarray-dataarray-sel) and slice to select a region of data\n",
    "- use [.mean](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.mean.html#xarray.DataArray.mean) to calculate the spatial mean\n",
    "- use [.plot](http://xarray.pydata.org/en/stable/generated/xarray.DataArray.plot.html#xarray-dataarray-plot) to plot a timeseries of the mean wind speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "subset = ds_ccmp.sel(latitude=slice(30,50),longitude=slice(200,210))\n",
    "\n",
    "ts = subset.mean({'latitude','longitude'},keep_attrs=True)\n",
    "\n",
    "ts.wspd.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a Hovmoller diagram\n",
    "\n",
    "- [Hovmoller diagrams](https://en.wikipedia.org/wiki/Hovm%C3%B6ller_diagram) are common plots in oceanography and meteorology to look at propogation of waves. They often have latitude or longitude along the x-axis and time along the y-axis so you can look at how things propogate in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "ds_ccmp.sel(latitude=0.125,longitude=slice(120,275)).wspd.plot(vmin=0,vmax=12,cmap='magma')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In Feb 2020 a GRL [paper](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2020GL091430) came out connecting 3 closely occuring Typhoons near Korea to the California wildfires\n",
    "  \n",
    "\"Strong winds that accentuated a fire outbreak in the western United States in early September of 2020 resulted from an atmospheric wave train that spanned the Pacific Ocean. Days before the atmospheric waves developed in the United States, three western Pacific tropical cyclones (typhoons) underwent an extratropical transition over Korea within an unprecedentedly short span of 12 days. \"\n",
    "\n",
    "Using ERA5 winds, Figure 1 showed zonal winds averaged over a box located over NCal/Oregon. On 9/8 and again on 9/24 the zonal winds were strongly negative (Diablo winds) and both events were associated with major fires.\n",
    "\n",
    "Below, let's do the same figure with CCMP winds - Note that CCMP winds end 10/4/2020, so the data will end at that point.\n",
    "\n",
    "- First create a timeseries of the 'uwnd' or east-west wind vector for the same region used in the paper above.\n",
    "- Next plot all the data in grey and then add 2020 data in red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = ds_ccmp.uwnd.sel(longitude=slice(236,239),latitude=slice(41,46)).mean({'latitude','longitude'}).load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot the data\n",
    "fig, ax = plt.subplots(figsize=(12, 6))\n",
    "for lyr in range(1979,2021):\n",
    "    yr = ts.sel(time=slice(str(lyr)+'-08-01',str(lyr)+'-10-29'))\n",
    "    xx=yr.time.dt.dayofyear+yr.time.dt.hour/24\n",
    "    ax.plot(xx,yr,color='silver')\n",
    "    if lyr==2020:\n",
    "        ax.plot(xx,yr,color='red')\n",
    "date_form = DateFormatter(\"%m-%d\")\n",
    "ax.xaxis.set_major_formatter(date_form)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare the satellite wind speeds to buoy (in situ) data\n",
    "Data from NDBC buoys, which measure wind speed are [here](https://dods.ndbc.noaa.gov/) and can be read via an THREDDS server.\n",
    "- read in NDBC buoy data\n",
    "- find closest CCMP data and linearly interpolate to buoy location\n",
    "- examine a timeseries, caluclate bias and STD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#read in NDBC buoy data\n",
    "url='https://dods.ndbc.noaa.gov/thredds/dodsC/data/cwind/51003/51003.ncml'\n",
    "\n",
    "with xr.open_dataset(url) as ds:\n",
    "    # The longitude is on a -180 to 180, CCMP is 0-360, so make sure to convert\n",
    "    ds.coords['longitude'] = np.mod(ds['longitude'], 360)\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#find the time period that both data exist eg. start and stop\n",
    "time_start = np.max([ds.time[0].data,ds_ccmp.time[0].data])\n",
    "time_stop = np.min([ds.time[-1].data,ds_ccmp.time[-1].data])\n",
    "\n",
    "#cut off satellite data to joint period\n",
    "# then linearly interpolate the data to the buoy location\n",
    "ccmp_buoy = ds_ccmp.sel(time=slice(time_start,time_stop)).interp(latitude=ds.latitude,longitude=ds.longitude,method='linear')\n",
    "\n",
    "#cut off buoy data to joint period\n",
    "ds = ds.sel(time=slice(time_start,time_stop))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# go from 30 min to 6-hourly sampling, resample with mean for +- 3 hours centered on timestamp\n",
    "#ds = ds.resample(time='6H',loffset='180min',base=3).mean()  \n",
    "# OMG!!! resample is sooo slow.  doing this cludge instead, same to within 6 significant digits\n",
    "\n",
    "ds = ds.rolling(time=36,center=True).mean() \n",
    "ds_buoy = ds.interp(time=ccmp_buoy.time,method='nearest')\n",
    "\n",
    "ds_buoy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#plot the data\n",
    "ccmp_buoy.wspd.plot(color='r')\n",
    "ds_buoy.wind_spd.plot(color='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# change the chunking in time for the satellite data\n",
    "ccmp_buoy = ccmp_buoy.chunk({\"time\": 24797}).compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# let's remove the seasonal cycle so we can compare the data more clearly\n",
    "\n",
    "#for ccmp data create the climatology and anomaly\n",
    "ts_ccmp_climatology = ccmp_buoy.groupby('time.dayofyear').mean('time',keep_attrs=True,skipna=True)\n",
    "ts_ccmp_anomaly = (ccmp_buoy.groupby('time.dayofyear')-ts_ccmp_climatology).compute()\n",
    "\n",
    "#for buoy data create the climatology and anomaly\n",
    "ts_buoy_climatology = ds_buoy.groupby('time.dayofyear').mean('time',keep_attrs=True,skipna=True)\n",
    "ts_buoy_anomaly = (ds_buoy.groupby('time.dayofyear')-ts_buoy_climatology).compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#plot the data\n",
    "(ts_ccmp_anomaly.wspd-ts_buoy_anomaly.wind_spd).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#calculate the mean difference and standard deviation\n",
    "\n",
    "print('buoy minus satellite wind speeds')\n",
    "rdif = (ts_buoy_anomaly.wind_spd-ts_ccmp_anomaly.wspd)\n",
    "print('mean:',rdif.mean().data)\n",
    "print('std:',rdif.std().data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Masking satellite data and looking at frequency of low/high wind speeds.\n",
    "- Make a land/ocean/ice mask to show where there is actually data\n",
    "\n",
    "#### Three different ways to mask the data\n",
    "1. A daily mask that removes data with sea ice and land\n",
    "- sum over time for nobs (number of observations) variable\n",
    "- average over a month so that land and monthly sea ice are masked out\n",
    "2. A mask that removes all data that over land or where there is 'permanent' sea ice\n",
    "- find when nobs is > 0\n",
    "3. A climatology mask that removes all data that over land or where there has ever been sea ice\n",
    "- sum over time for nobs (number of observations) variable\n",
    "- average over a month so that land and monthly sea ice are masked out\n",
    "\n",
    "#### Apply the mask \n",
    "- over land, CCMP is ERA5 data\n",
    "- for many ocean applications a land / sea ice mask is needed\n",
    "- below are some different mask options that use the CCMP data to generate a mask\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mask_data(ds,type):\n",
    "    if type=='daily': #daily mask removes sea ice and land\n",
    "        mask_obs = ds.nobs.rolling(time=180,center=True).max('time')  #4 per day 30 days = 180 rolling window\n",
    "        cutoff = 0\n",
    "    if type=='land':  # land mask only (includes data over sea ice)\n",
    "        mask_obs = ds.nobs.sum({'time'},keep_attrs=True)  #this will give you a LAND mask\n",
    "        cutoff = 0\n",
    "    if type=='climatology':  #climatology mask removes max sea ice extent and land\n",
    "        test = ds.nobs.coarsen(time=16, boundary=\"trim\").max() #subsample to every 4 days\n",
    "        test = test/test #normalize\n",
    "        mask_obs = test.rolling(time=10,center=True).max('time')\n",
    "        mask_obs = mask_obs.sum({'time'},keep_attrs=True).compute()\n",
    "        cutoff = mask_obs.sel(longitude=slice(220,230),latitude=slice(10,20)).min() #pick area in middle pacific to get #obs\n",
    "        cutoff = cutoff-cutoff*.2 #just because there are some areas with a little less data eg. itcz\n",
    "    dy_mask = (mask_obs>cutoff).compute()  #computing the mask speeds up subsequent operations\n",
    "    masked = ds.where(dy_mask)\n",
    "    return masked,dy_mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For this tutorial we will use the climatology mask\n",
    "\n",
    "- This code will take ~10 min to run through the entire dataset to create the climatology mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "subset = ds_ccmp.sel(time='2018')\n",
    "\n",
    "masked,mask_obs = mask_data(subset,'climatology')\n",
    "\n",
    "mask_obs.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(18,6))\n",
    "\n",
    "masked.wspd[100,:,:].plot(ax=ax[0])\n",
    "\n",
    "masked.wspd[-100,:,:].plot(ax=ax[1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When are there low or high wind speeds at each location?\n",
    "\n",
    "- calculate percentage of low/high winds at each location\n",
    "\n",
    "- show a spatial map showing a climatology of roughness.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from xhistogram.xarray import histogram\n",
    "\n",
    "# wspd_bin = 0 gives fraction of winds <= 2 m/s where strong diurnal warming of the ocean surface may be present\n",
    "\n",
    "wspd_levels = np.array([0,2,1000]) \n",
    "\n",
    "wspd_hist = histogram(subset.wspd, \n",
    "                      bins=[wspd_levels], \n",
    "                      dim=['time'], \n",
    "                      density=False)\n",
    "\n",
    "wspd_hist = wspd_hist/wspd_hist.sum('wspd_bin')\n",
    "\n",
    "wspd_hist.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from xhistogram.xarray import histogram\n",
    "\n",
    "# wspd_bin = 1 gives fraction of winds > 10 m/s where wind speed accuracy decreases\n",
    "\n",
    "wspd_levels = np.array([0,10,1000]) \n",
    "wspd_hist2 = histogram(subset.wspd, \n",
    "                      bins=[wspd_levels], \n",
    "                      dim=['time'], \n",
    "                      density=False)\n",
    "\n",
    "wspd_hist2 = wspd_hist2/wspd_hist2.sum('wspd_bin')\n",
    "\n",
    "wspd_hist2.load()\n",
    "\n",
    "whist = xr.concat([wspd_hist.isel(wspd_bin=0),wspd_hist2.isel(wspd_bin=1)],dim='wspd_bin')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "#plot the data\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (8.0,4.0)\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "fg = whist.plot(vmin=0, vmax=1,\n",
    "    col=\"wspd_bin\",\n",
    "    transform=ccrs.PlateCarree(),  # remember to provide this!\n",
    "    subplot_kws={\n",
    "        \"projection\": ccrs.PlateCarree()\n",
    "    },\n",
    "    cbar_kwargs={\"label\":'Fraction of obs',\"orientation\": \"horizontal\", \"shrink\": 0.8, \"aspect\": 40},\n",
    "    robust=True,\n",
    ")\n",
    "tstr = ['< 2 m/s','> 10 m/s']\n",
    "for i, ax in enumerate(fg.axes.flat):\n",
    "    ax.set_title(tstr[i]) \n",
    "fg.map(lambda: plt.gca().coastlines())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate using only data over the ocean\n",
    "- use the mask to remove the land data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "masked = whist.where(mask_obs>0)\n",
    "\n",
    "fg = masked.plot(vmin=0, vmax=1,\n",
    "    col=\"wspd_bin\",\n",
    "    transform=ccrs.PlateCarree(),  # remember to provide this!\n",
    "    subplot_kws={\n",
    "        \"projection\": ccrs.PlateCarree()\n",
    "    },\n",
    "    cbar_kwargs={\"label\":'Fraction of obs',\"orientation\": \"horizontal\", \"shrink\": 0.8, \"aspect\": 40},\n",
    "    robust=True,\n",
    ")\n",
    "tstr = ['< 2 m/s','> 10 m/s']\n",
    "for i, ax in enumerate(fg.axes.flat):\n",
    "    ax.set_title(tstr[i]) \n",
    "fg.map(lambda: plt.gca().coastlines())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client.close()"
   ]
  }
 ],
 "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}