{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# # YOPP Forecast\n", "\n", "# - Loads in all daily forecasts of sea ice extent\n", "# - Regrids to polar stereographic,\n", "# - Saves to netcdf files grouped by year" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "'''\n", "\n", "This code is part of the SIPN2 project focused on improving sub-seasonal to seasonal predictions of Arctic Sea Ice. \n", "If you use this code for a publication or presentation, please cite the reference in the README.md on the\n", "main page (https://github.com/NicWayand/ESIO). \n", "\n", "Questions or comments should be addressed to nicway@uw.edu\n", "\n", "Copyright (c) 2018 Nic Wayand\n", "\n", "GNU General Public License v3.0\n", "\n", "\n", "'''\n", "\n", "# Standard Imports\n", "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload\n", "import matplotlib\n", "import scipy\n", "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", "from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import xesmf as xe\n", "import os\n", "import glob\n", "import seaborn as sns\n", "import warnings\n", "warnings.simplefilter(action='ignore', category=FutureWarning)\n", "\n", "# ESIO Imports\n", "\n", "from esio import EsioData as ed\n", "from esio import import_data" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# General plotting settings\n", "sns.set_style('whitegrid')\n", "sns.set_context(\"talk\", font_scale=1.5, rc={\"lines.linewidth\": 2.5})" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "E = ed.EsioData.load()\n", "# Directories\n", "model='yopp'\n", "runType='forecast'\n", "updateall = False\n", "data_dir = E.model[model][runType]['native']\n", "data_out = E.model[model][runType]['sipn_nc']\n", "model_grid_file = E.model[model]['grid']\n", "stero_grid_file = E.obs['NSIDC_0051']['grid']" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "obs_grid = import_data.load_grid_info(stero_grid_file, model='NSIDC')\n", "# Ensure latitude is within bounds (-90 to 90)\n", "# Have to do this because grid file has 90.000001\n", "obs_grid['lat_b'] = obs_grid.lat_b.where(obs_grid.lat_b < 90, other = 90)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Regridding Options\n", "method='nearest_s2d' # ['bilinear', 'conservative', 'nearest_s2d', 'nearest_d2s', 'patch']" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## TODO\n", "# - Get mask\n", "# - Get lat lon bounds " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "all_files = glob.glob(os.path.join(data_dir, 'yopp*ci*.grib'))\n", "print(\"Found \",len(all_files),\" files.\")\n", "if updateall:\n", " print(\"Updating all files...\")\n", "else:\n", " print(\"Only updating new files\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "weights_flag = False # Flag to set up weights have been created\n", "\n", "cvar = 'sic'\n", "\n", "# Load land/sea mask file\n", "if os.path.basename(model_grid_file)!='MISSING':\n", " ds_mask = xr.open_dataset(model_grid_file)\n", "else:\n", " ds_mask = None\n", "\n", "for cf in all_files:\n", " # Check if already imported and skip (unless updateall flag is True)\n", " f_out = os.path.join(data_out, os.path.basename(cf).split('.')[0]+'_Stereo.nc') # netcdf file out \n", " if not updateall:\n", " if os.path.isfile(f_out):\n", " print(\"Skipping \", os.path.basename(cf), \" already imported.\")\n", " continue # Skip, file already imported\n", "\n", " ds = xr.open_dataset(cf, engine='pynio')\n", "\n", " # Rename variables per esipn guidelines\n", " ds.rename({'CI_GDS4_SFC':'sic', 'g4_lat_2':'lat', 'g4_lon_3':'lon', 'initial_time0_hours':'init_time',\n", " 'forecast_time1':'fore_time'}, inplace=True);\n", " \n", " # Apply masks (if available)\n", " if ds_mask:\n", " print('found mask')\n", " # land_mask is the fraction of native grid cell that is land\n", " # (1-land_mask) is fraction ocean\n", " # Multiply sic by fraction ocean to get actual native grid cell sic\n", " # Also mask land out where land_mask==1\n", " ds[cvar] = ds[cvar] * (1 - ds_mask.land_mask.where(ds_mask.land_mask<0.5)) # Using 50% threshold here\n", " \n", "\n", "# ds.coords['nj'] = model_grid.nj\n", "# ds.coords['ni'] = model_grid.ni\n", "# ds.coords['lat'] = model_grid.lat\n", "# ds.coords['lon'] = model_grid.lon\n", "# ds.coords['lat_b'] = model_grid.lat_b\n", "# ds.coords['lon_b'] = model_grid.lon_b\n", "# ds.coords['imask'] = model_grid.imask\n", " \n", "# # Set sic below 0 to 0\n", "# if X.sic.min().values < 0:\n", "# print(\"Some negative SIC \"+str(X.sic.min().values)+\" found in input PIOMAS, setting to 0\")\n", "# X = X.where(X>=0, other=0)\n", " \n", "# # Apply model mask\n", "# X = X.where(X.imask)\n", " \n", " # Calculate regridding matrix\n", " regridder = xe.Regridder(ds, obs_grid, method, periodic=False, reuse_weights=weights_flag)\n", " weights_flag = True # Set true for following loops\n", " \n", " # Add NaNs to empty rows of matrix (forces any target cell with ANY source cells containing NaN to be NaN)\n", " if method=='conservative':\n", " regridder = import_data.add_matrix_NaNs(regridder)\n", " \n", " # Regrid variable\n", " var_out = regridder(ds[cvar])\n", " \n", " # Expand dims\n", " var_out = import_data.expand_to_sipn_dims(var_out)\n", "\n", " # # Save regridded to netcdf file\n", " \n", " var_out.to_netcdf(f_out)\n", " var_out = None # Memory clean up\n", " print('Saved ', f_out)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Clean up\n", "if weights_flag:\n", " regridder.clean_weight_file() # clean-up" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "# sic_all = xr.open_dataset(f_out)\n", "\n", "\n", "\n", "\n", "\n", "# # Set up plotting info\n", "# cmap_sic = matplotlib.colors.ListedColormap(sns.color_palette(\"Blues\", 10))\n", "# cmap_sic.set_bad(color = 'red')\n", "\n", "# # Plot original projection\n", "# plt.figure(figsize=(20,10))\n", "# ax1 = plt.axes(projection=ccrs.PlateCarree())\n", "# ds_p = ds.sic.isel(init_time=1).isel(fore_time=79)\n", "# ds_p.plot.pcolormesh(ax=ax1, x='lon', y='lat', \n", "# vmin=0, vmax=1,\n", "# cmap=matplotlib.colors.ListedColormap(sns.color_palette(\"Blues\", 10)),\n", "# transform=ccrs.PlateCarree());\n", "# ax1.set_extent([-180, 180, -90, 90], crs=ccrs.PlateCarree())\n", "# gl = ax1.gridlines(crs=ccrs.PlateCarree(), linestyle='-')\n", "# gl.xlabels_bottom = True\n", "# gl.ylabels_left = True\n", "# gl.xformatter = LONGITUDE_FORMATTER\n", "# gl.yformatter = LATITUDE_FORMATTER\n", "# ax1.coastlines(linewidth=0.75, color='black', resolution='50m');\n", "\n", "\n", "# # Plot SIC on target projection\n", "# (f, ax1) = ice_plot.polar_axis()\n", "# f.set_size_inches((10,10))\n", "# ds_p.plot.pcolormesh(ax=ax1, x='lon', y='lat', \n", "# transform=ccrs.PlateCarree(),\n", "# cmap=cmap_sic)\n", "# ax1.set_title('Original Grid')\n", "\n", "\n", "# # Plot SIC on target projection\n", "# (f, ax1) = ice_plot.polar_axis()\n", "# f.set_size_inches((10,10))\n", "# ds_p2 = sic_all.sic.isel(init_time=1).isel(fore_time=79).isel(ensemble=0)\n", "# ds_p2.plot.pcolormesh(ax=ax1, x='lon', y='lat', \n", "# transform=ccrs.PlateCarree(),\n", "# cmap=cmap_sic)\n", "# ax1.set_title('Target Grid')\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "py36 pynio", "language": "python", "name": "test_nio" }, "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }