notebooks/Import_NSIDC_Obs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Importing  NSIDC_0081 ...\n",
      "Found  1  new files to import...\n",
      "Finished  NSIDC_0081\n",
      "\n",
      "Importing  NSIDC_0051 ...\n",
      "Found  0  new files to import...\n",
      "Finished  NSIDC_0051\n",
      "\n",
      "Importing  NSIDC_0079 ...\n",
      "Found  0  new files to import...\n",
      "Finished  NSIDC_0079\n",
      "\n"
     ]
    }
   ],
   "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",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import os\n",
    "import xarray as xr\n",
    "import glob\n",
    "# import loadobservations as lo\n",
    "from esio import import_data\n",
    "from esio import metrics\n",
    "from esio import EsioData as ed\n",
    "\n",
    "# Dirs\n",
    "E = ed.EsioData.load()\n",
    "data_dir = E.obs_dir\n",
    "\n",
    "# Flags\n",
    "UpdateAll = False\n",
    "\n",
    "# Load in regional data\n",
    "# Note minor -0.000004 degree differences in latitude\n",
    "ds_region = xr.open_dataset(os.path.join(E.grid_dir, 'sio_2016_mask_Update.nc'))\n",
    "\n",
    "# Products to import\n",
    "product_list = ['NSIDC_0081', 'NSIDC_0051', 'NSIDC_0079']\n",
    "\n",
    "# Version numbers (if any)\n",
    "# If multiple version number files exist, it only imports the one specified below\n",
    "ver_nums = {'NSIDC_0079':'v3.1','NSIDC_0081':'nrt','NSIDC_0051':'v1.1'}\n",
    "\n",
    "ds_lat_lon = import_data.get_stero_N_grid(grid_dir=E.grid_dir)\n",
    "\n",
    "# Loop through each product\n",
    "for c_product in product_list:\n",
    "    print('Importing ', c_product, '...')\n",
    "\n",
    "    # Find new files that haven't been imported yet\n",
    "    native_dir = os.path.join(data_dir, c_product, 'native')\n",
    "    os.chdir(native_dir)\n",
    "    native_files = sorted(glob.glob('*'+ver_nums[c_product]+'*.bin'))\n",
    "    nc_dir = os.path.join(data_dir, c_product, 'sipn_nc')\n",
    "    os.chdir(nc_dir)\n",
    "    nc_files = sorted(glob.glob('*.nc'))\n",
    "    if UpdateAll:\n",
    "        new_files = [x.split('.b')[0] for x in native_files]\n",
    "        print('Updating all ', len(native_files), ' files...')\n",
    "    else:\n",
    "        new_files = np.setdiff1d([x.split('.b')[0] for x in native_files], \n",
    "                                 [x.split('.n')[0] for x in nc_files]) # just get file name and compare\n",
    "        print('Found ', len(new_files), ' new files to import...')\n",
    "\n",
    "    # Loop through each file\n",
    "    for nf in new_files:\n",
    "        \n",
    "        # Load in \n",
    "        ds_sic = import_data.load_1_NSIDC(filein=os.path.join(native_dir, nf+'.bin'), product=c_product)\n",
    "\n",
    "        # Add lat and lon dimensions\n",
    "        ds_sic.coords['lat'] = ds_lat_lon.lat\n",
    "        ds_sic.coords['lon'] = ds_lat_lon.lon\n",
    "\n",
    "        # Stereo projected units (m)\n",
    "        dx = dy = 25000 \n",
    "        xm = np.arange(-3850000, +3750000, +dx)\n",
    "        ym = np.arange(+5850000, -5350000, -dy)\n",
    "        ds_sic.coords['xm'] = xr.DataArray(xm, dims=('x'))\n",
    "        ds_sic.coords['ym'] = xr.DataArray(ym, dims=('y'))    \n",
    "        \n",
    "        # Calculate extent and area\n",
    "#         ds_sic['extent'] = ((ds_sic.sic>=0.15).astype('int') * ds_region.area).sum(dim='x').sum(dim='y')/(10**6)\n",
    "        ds_sic['extent'] = metrics.calc_extent(ds_sic.sic, ds_region, fill_pole_hole=True)\n",
    "#         ds_sic['extent'] = ds_sic['extent'] + (ds_sic.hole_mask.astype('int') * ds_region.area).sum(dim='x').sum(dim='y')/(10**6) # Add hole\n",
    "        ds_sic['area'] = (ds_sic.sic * ds_region.area).sum(dim='x').sum(dim='y')/(10**6) # No pole hole\n",
    "    \n",
    "        # Save to netcdf file\n",
    "        ds_sic.to_netcdf(os.path.join(nc_dir, nf.split('.b')[0]+'.nc'))\n",
    "        ds_sic = None\n",
    "    \n",
    "#     # Calculate extent and area (saved to separte file)\n",
    "#     if len(new_files) > 0 : # There were some new files\n",
    "#         print('Calculating extent and area...')\n",
    "#         ds_all = xr.open_mfdataset(os.path.join(nc_dir,'*.nc'), concat_dim='time', \n",
    "#                                    autoclose=True, compat='no_conflicts',\n",
    "#                                    data_vars=['sic'])    \n",
    "#         print('Loaded in all files...')\n",
    "#         ds_all['extent'] = ((ds_all.sic>=0.15).astype('int') * ds_region.area).sum(dim='x').sum(dim='y')/(10**6)\n",
    "#         ds_all['extent'] = ds_all['extent'] + (ds_all.hole_mask.astype('int') * ds_region.area).sum(dim='x').sum(dim='y')/(10**6) # Add hole\n",
    "#         ds_all['area'] = (ds_all.sic * ds_region.area).sum(dim='x').sum(dim='y')/(10**6) # No pole hole\n",
    "#         ds_all = ds_all[['extent','area']]\n",
    "#         # Create new dir to store agg file\n",
    "#         if not os.path.exists(os.path.join(data_dir, c_product, 'agg_nc')):\n",
    "#             os.makedirs(os.path.join(data_dir, c_product, 'agg_nc'))\n",
    "#         ds_all.to_netcdf(os.path.join(data_dir, c_product, 'agg_nc', 'panArctic.nc'))\n",
    "        \n",
    "    # For each Product\n",
    "    print(\"Finished \", c_product)\n",
    "    print(\"\")"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.6.4 esio",
   "language": "python",
   "name": "esio"
  },
  "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
}