code/scripts/h5_to_np_icesat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python 3 script to extract geophysical data from ICESAT HDF5 format files into python numpy-based arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 0 Load libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "import dask\n",
    "import dask.dataframe as dd\n",
    "from dask.diagnostics import ProgressBar\n",
    "import glob\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import tables\n",
    "import xarray as xr\n",
    "\n",
    "print(os.getcwd())\n",
    "print(_dh[0])\n",
    "os.chdir(os.environ['HOME']+\"/data/icesat/GLAH12.034/2003.02.20\")\n",
    "print(os.getcwd())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1 Load Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declare functions that stores the fields we want to retrieve from ICESAT (e.g. xyzt data) and then extract those fields from the HDF5 file into Python formats like Numpy/Pandas/Xarray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_h5_keyDict(h5file, datagroup=\"Data_40HZ\", useAll=False):\n",
    "    '''\n",
    "    Function to pick the fields we want from a HDF5 file.\n",
    "    Currently hardcoded for ICESAT GLAH12 HDF5 files.\n",
    "\n",
    "    Arguments:\n",
    "    h5file -- input data, of type HDF5 https://support.hdfgroup.org/HDF5/whatishdf5.html\n",
    "    datagroup -- the HDF5 data group to look in for ICESAT see https://nsidc.org/data/glas/data-dictionary-glah12\n",
    "    useAll -- Set to True if you want to retrieve all fields in the datagroup, else just get stardard xyzt data\n",
    "\n",
    "    Returns:\n",
    "    fields -- python dictionary with keys (e.g. x, y, z, t) mapped to values (fieldnames within the HDF5 file)\n",
    "    '''\n",
    "    h5 = h5py.File(h5file, \"r\")\n",
    "    assert(isinstance(h5, (h5py._hl.files.File)))\n",
    "\n",
    "    #GLAH12 Product Data Dictionary https://nsidc.org/data/glas/data-dictionary-glah12\n",
    "    #[k for k in h5.keys()]\n",
    "\n",
    "    #%% Retrieve data stored in the hdf5 file using known keys.\n",
    "    [g for g in h5['{0}'.format(datagroup)]]\n",
    "\n",
    "    fields = collections.OrderedDict()\n",
    "    [v.name for v in h5[datagroup].values()]\n",
    "\n",
    "    #Standard parameters\n",
    "    if useAll == False:\n",
    "        fields['i'] = datagroup+'/Time/i_rec_ndx'                 #GLAS Record Index (i)\n",
    "        fields['x'] = datagroup+'/Geolocation/d_lon'              #Longitude         (x)\n",
    "        fields['y'] = datagroup+'/Geolocation/d_lat'              #Latitude          (y)\n",
    "        if datagroup == 'Data_40HZ':\n",
    "            fields['z'] = datagroup+'/Elevation_Surfaces/d_elev'  #Elevation         (z)\n",
    "            fields['t'] = datagroup+'/Time/d_UTCTime_40'          #Timestamp         (t)\n",
    "        elif datagroup == 'Data_1HZ':\n",
    "            fields['k'] = datagroup+'/Geolocation/i_track'        #Track number      (k)\n",
    "            fields['t'] = datagroup+'/Time/d_UTCTime_1'           #Timestamp         (t)\n",
    "\n",
    "    #All other useful-ish parameters\n",
    "    if useAll == True:\n",
    "        def func(name, obj):\n",
    "            if isinstance(obj, h5py.Dataset):\n",
    "                if obj.ndim == 1:\n",
    "                    fields[name] = datagroup+'/'+name\n",
    "                else:\n",
    "                    print(\"Warn: {0} is not one-dimensional, ignoring...\".format(name))\n",
    "        h5[datagroup].visititems(func)\n",
    "\n",
    "    assert(isinstance(fields, dict))\n",
    "    return fields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def h5_to_pydata(h5file, h5fields):\n",
    "    '''\n",
    "    Function to load data from one HDF5 file to standard Python based formats.\n",
    "    Basically creates an n-dimensional array to store all the data. Do this by stacking each feature (e.g. coordinates/time/other param) using np.stack.\n",
    "\n",
    "    Arguments:\n",
    "    h5file -- input data, of type HDF5 https://support.hdfgroup.org/HDF5/whatishdf5.html\n",
    "    h5fields -- python dictionary containing keys (e.g. x, y, z, t) mapped to values (fieldnames within the HDF5 file)\n",
    "\n",
    "    Returns:\n",
    "    npData  -- numpy.array of shape (n, m) where n is number of fields and m is number of datapoints\n",
    "    pdData -- pandas.DataFrame version of above numpy array, with the time (t) field/column expressed in standard python datetime format\n",
    "    xrData -- xarray.Dataset version of the above pandas.DataFrame which has dimensions (m) and data variables (n)\n",
    "    '''\n",
    "    h5 = h5py.File(h5file, \"r\")\n",
    "    assert(isinstance(h5, (h5py._hl.files.File)))\n",
    "    #h5.name\n",
    "    #h5.libver\n",
    "    #h5.driver\n",
    "\n",
    "    #%% calculate m (number of individual datapoints) so we can do reshapes and assertion checks\n",
    "    dataListShape = [h5[h5fields[key]].shape for key in h5fields.keys()]\n",
    "    assert(np.median(dataListShape) == np.max(dataListShape))  #stupid way to get 'm' which is the no. of individual datapoints\n",
    "    m = np.max(dataListShape) #take m as the largest length\n",
    "\n",
    "    #%% numpy\n",
    "    npData = np.hstack((h5[h5fields[key]][:].reshape(-1,1) for key in h5fields.keys() if h5[h5fields[key]].shape == m)).T\n",
    "    assert(npData.shape == (len(h5fields), m))  #check that final numpy array has shape (n, m) where n is no. of features and m is no. of datapoints e.g. (4, 20000)\n",
    "    npData.shape\n",
    "    npData.ndim\n",
    "    npData.T.ndim\n",
    "\n",
    "    #%% pandas\n",
    "    assert(isinstance(npData, np.ndarray))\n",
    "    pdData = pd.DataFrame(npData.T, columns=h5fields.keys())\n",
    "    pdData['t'] = pd.to_datetime(pdData['t'], unit='s', origin=pd.Timestamp('2000-01-01'), infer_datetime_format=True)  #convert time data into standard python datetime format\n",
    "    assert(isinstance(pdData['t'][0], pd.Timestamp))\n",
    "\n",
    "    #%% xarray\n",
    "    assert(isinstance(pdData, pd.DataFrame))\n",
    "    xrData = pdData.to_xarray()\n",
    "    xrData\n",
    "    assert(isinstance(xrData, xr.Dataset))\n",
    "\n",
    "    return npData, pdData, xrData"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h5fields40hz = init_h5_keyDict(\"GLAH12_634_1102_001_0071_0_01_0001.H5\", datagroup=\"Data_40HZ\", useAll=False)\n",
    "h5fields1hz = init_h5_keyDict(\"GLAH12_634_1102_001_0071_0_01_0001.H5\", datagroup=\"Data_1HZ\", useAll=False)\n",
    "\n",
    "npData, pdData, xrData = h5_to_pydata(\"GLAH12_634_1102_001_0071_0_01_0001.H5\", h5fields40hz)\n",
    "df40 = pdData.loc[:,['x','y','z','i']].loc[lambda df: df.y < 0]  #filter for Antarctica only (South of Equator)\n",
    "\n",
    "npData, pdData, xrData = h5_to_pydata(\"GLAH12_634_1102_001_0071_0_01_0001.H5\", h5fields1hz)\n",
    "df1 = pdData.loc[:,['x','y','k','i']].loc[lambda df: df.y < 0]  #filter for Antarctica only (South of Equator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Export the data to csv format if we so wish to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df40\n",
    "#pdData.to_csv(os.environ['HOME']+\"/code/scripts/pdData.csv\")   #export Greenland and Antarctic data\n",
    "#df.to_csv(os.environ['HOME']+\"/code/scripts/pdData.csv\")       #export Antarctic data (South of Equator) only"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Dask\n",
    "#See also VITables, a GUI for PyTables https://github.com/uvemas/ViTables\n",
    "print(dask.__version__)\n",
    "p = ProgressBar()  #Real-time feedback on dask processes\n",
    "p.register()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Info inside the HDF5 file\n",
    "#h5file = os.path.join(os.environ['HOME']+\"/data/icesat/GLAH12.034/2003.02.20/GLAH12_634_1102_001_0071_0_01_0001_test.H5\")\n",
    "#[k for k in h5py.File(h5file, \"r\").keys()]\n",
    "#list(h5fields40hz.values())[0]\n",
    "\n",
    "#%% Debugging https://github.com/pandas-dev/pandas/issues/17661\n",
    "#store = pd.HDFStore(h5file, mode='r+')\n",
    "#store.select('Data_40HZ/Geolocation/d_lon')\n",
    "#store.__contains__(\"Data_40HZ/Geolocation\")\n",
    "#store.get_node(\"Data_40HZ/Elevation_Surfaces\")\n",
    "#store.get_storer(\"Data_40HZ\")\n",
    "#store.groups()\n",
    "#store.items\n",
    "#store.close()\n",
    "\n",
    "#tables.is_hdf5_file(h5file)\n",
    "#lala = tables.open_file(h5file, mode='r+', root_uep='/Data_40HZ/Geolocation')\n",
    "#lala.del_node_attr('/d_lon', 'DIMENSION_LIST')\n",
    "#lala.close()\n",
    "\n",
    "#pd.read_hdf(h5file, key='/Data_40HZ/Geolocation/d_lon')\n",
    "#dd.read_hdf(h5file, key='Data_40HZ/Geolocation/d_lon')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Get dask.dataframe(s) for fields [i, x, y, z, k, t] from Data_40HZ and Data_1HZ\n",
    "bugFixed = False\n",
    "if bugFixed == True:\n",
    "    #ideal command to run once ICESAT dask.read_hdf bug is fixed, upstream problem with PyTables https://github.com/PyTables/PyTables/issues/647\n",
    "    df40 = dd.read_hdf(os.environ['HOME']+\"/data/icesat/GLAH12.034/**/*.H5', key='/Data_40HZ\")\n",
    "    df1 = dd.read_hdf(os.environ['HOME']+\"/data/icesat/GLAH12.034/**/*.H5', key='/Data_1HZ\")\n",
    "else:\n",
    "    hpyPath = os.environ['HOME']+\"/data/icesat/GLAHPY12.034\"\n",
    "    os.makedirs(hpyPath, mode=0o777, exist_ok=True)\n",
    "    def pdData_to_hdf(h5f):\n",
    "        outFile = hpyPath+\"/\"+h5f.split('/')[-1]\n",
    "        if not os.path.exists(outFile):\n",
    "            print(h5f.split('/')[-1])\n",
    "            pd40 = h5_to_pydata(h5f, h5fields40hz)[1]\n",
    "            pd40.to_hdf(outFile, key=\"/Data_40HZ\", format='table', mode='w')\n",
    "            pd1 = h5_to_pydata(h5f, h5fields1hz)[1]\n",
    "            pd1.to_hdf(outFile, key=\"/Data_1HZ\", format='table', mode='a')\n",
    "    if len(glob.glob(hpyPath+'/*.H5')) != 637:\n",
    "        [pdData_to_hdf(h) for h in glob.iglob(os.environ['HOME']+\"/data/icesat/GLAH12.034/**/*.H5\")]  #convert data from raw NSIDC supplied HDF5 into PyTables compatible format\n",
    "    subsetLen = 20\n",
    "    subsetFiles = glob.glob(hpyPath+'/*.H5')[:subsetLen]\n",
    "    df40 = dd.read_hdf(subsetFiles, key='/Data_40HZ')  #workaround command to load Data_40HZ data into dask, not dask.delayed so slow\n",
    "    df1 = dd.read_hdf(subsetFiles, key='/Data_1HZ')   #workaround command to load Data_40HZ data into dask, not dask.delayed so slow\n",
    "    #df40 = dask.delayed(dd.read_hdf)(hpyPath+'/*.H5', key='/Data_40HZ')\n",
    "    #df1 = dask.delayed(dd.read_hdf)(hpyPath+'/*.H5', key='/Data_1HZ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#%% Join the Data_40HZ and Data_1HZ data on i_rec_ndx(i) to get Track Number(k) on Data_40HZ table\n",
    "assert(list(h5fields40hz.keys()) == ['i', 'x', 'y', 'z', 't'])\n",
    "assert(list(h5fields1hz.keys()) == ['i', 'x', 'y', 'k', 't'])\n",
    "df_all = dask.delayed(df40.merge)(df1[['i', 'k']], on='i')   #Perform dask delayed parallel join\n",
    "os.chdir(os.environ['HOME']+\"/code/scripts\")                      #change directory so that mydask.png can be saved to the right directory when running dask.visualize()\n",
    "#df_all.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#%% Computationally intensive code if running on full ICESAT dataset!!\n",
    "df = df_all.compute()   #very computationally intensive!!\n",
    "assert(isinstance(df, dd.DataFrame))\n",
    "#df.to_csv(hpyPath+'/export*.csv')  #export the joined table into csv files\n",
    "\n",
    "#df_all['k'].unique().compute()   #compute unique values of 'k' where k is the ICESAT Track Number\n",
    "toFilter = True\n",
    "if toFilter == True:\n",
    "    trackSubset = ((df_all['k'] <= 72) | (df_all['k'] >= 42))  #which ICESAT tracks to use\n",
    "    boundSubset = ((df_all['y'] < 0) & (df_all['x'] >= 0))   #Geographical boundaries\n",
    "    sqlFilter = trackSubset & boundSubset\n",
    "    df = df_all[sqlFilter]\n",
    "elif toFilter == False:\n",
    "    df = df_all\n",
    "#df.visualize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not isinstance(df, dd.DataFrame):\n",
    "    df = df.compute()  #Computes lazy dataframe and makes it non-lazy\n",
    "    #df = df.persist()  #will persist dataframe in RAM\n",
    "print(df.__class__, df)\n",
    "assert(isinstance(df, dd.DataFrame))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2 Plot those datapoints!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%matplotlib notebook\n",
    "%matplotlib inline\n",
    "#!conda install -y -c conda-forge nodejs\n",
    "#!jupyter labextension install jupyterlab_bokeh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Holoviews + Datashader + Geoviews"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!conda install -y holoviews\n",
    "#!conda install -y -c bokeh datashader\n",
    "#!conda install -c ioam geoviews\n",
    "import holoviews as hv\n",
    "import datashader as ds\n",
    "import geoviews as gv\n",
    "print(hv.__version__)\n",
    "print(ds.__version__)\n",
    "print(gv.__version__)\n",
    "hv.extension('bokeh')\n",
    "hv.notebook_extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#holoviews misc imports\n",
    "from holoviews.streams import * #RangeXY\n",
    "from colorcet import cm\n",
    "#datashader misc imports\n",
    "from holoviews.operation.datashader import aggregate, datashade, dynspread, shade\n",
    "#Geoviews misc imports\n",
    "import geoviews.feature as gf\n",
    "from cartopy import crs\n",
    "#Datashader options\n",
    "dynspread.max_px=20\n",
    "dynspread.threshold=0.5\n",
    "shade.cmap=\"#30a2da\" # to match HV Bokeh default\n",
    "# See https://anaconda.org/jbednar/holoviews_datashader/notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(df.info(), df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#df = df.set_index('t')\n",
    "#Reproject points from EPSG 4326 (PlateCarree) to EPSG 3031\n",
    "points = gv.Points(df, kdims=[('x', '3031X'), ('y', '3031Y'), ('z', 'Altitude')], vdims=['i', 'k'], crs = crs.PlateCarree())\n",
    "projected_gv = gv.operation.project_points(points, projection=crs.epsg(3031))\n",
    "assert(isinstance(projected_gv, gv.element.geo.Points))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%opts QuadMesh [tools=['hover']] (alpha=0 hover_alpha=0.2)\n",
    "%%output size=300  #set output size, e.g. 200 = 2x the default output size\n",
    "hv.notebook_extension('bokeh')\n",
    "gv_options = {'bgcolor':'black', 'show_grid':True}\n",
    "\n",
    "hvmap = projected_gv\n",
    "dsmap = datashade(hvmap, x_sampling=1, y_sampling=1, cmap=cm['fire'])\n",
    "gvmap = dynspread(dsmap.opts(plot=gv_options))\n",
    "#gvmap * hv.util.Dynamic(aggregate(hvmap, width=5, height=5, streams=[PointerX]), operation=hv.QuadMesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#!pip install paramnb\n",
    "#!conda install -y -c ioam parambokeh\n",
    "import param\n",
    "import paramnb\n",
    "#import parambokeh\n",
    "\n",
    "minAlt = round(df['z'].min().compute(), -2)\n",
    "maxAlt = round(df['z'].max().compute(), -2)\n",
    "minTrk = int(df['k'].min().compute())\n",
    "maxTrk = int(df['k'].max().compute())\n",
    "\n",
    "\n",
    "class IceSatExplorer(hv.streams.Stream):\n",
    "    colormap   = param.ObjectSelector(default=cm[\"fire\"], objects=cm.values())\n",
    "    altitude   = param.Range(default=(minAlt, maxAlt), bounds=(minAlt, maxAlt), doc=\"\"\"Elevation of ICESAT laser point\"\"\")\n",
    "    trackrange = param.Range(default=(minTrk, maxTrk), bounds=(minTrk, maxTrk), doc=\"\"\"ICESAT Track subset\"\"\" )\n",
    "    timerange  = param.Range()\n",
    "    \n",
    "    def make_view(self, x_range=None, y_range=None, **kwargs):\n",
    "        #map_tiles = tiles.opts(style=dict(alpha=self.alpha), plot=options) \n",
    "\n",
    "        hvmap = projected_gv.select(z=self.altitude, k=self.trackrange)\n",
    "        dsmap = datashade(hvmap, x_sampling=1, y_sampling=1, cmap=self.colormap,\n",
    "                          dynamic=False, x_range=x_range, y_range=y_range)\n",
    "        gv_options = {'bgcolor':'black', 'show_grid':True}\n",
    "        gvmap = dynspread(dsmap.opts(plot=gv_options))\n",
    "        return gvmap #* hv.util.Dynamic(aggregate(hvmap, width=5, height=5, streams=[PointerX]), operation=hv.QuadMesh)\n",
    "        \n",
    "        #points = hv.Points(df, kdims=[self.plot+'_x', self.plot+'_y'], vdims=['passenger_count'])\n",
    "        #selected = points.select(passenger_count=self.passengers)\n",
    "        #taxi_trips = datashade(selected, x_sampling=1, y_sampling=1, cmap=self.colormap,\n",
    "        #                       dynamic=False, x_range=x_range, y_range=y_range,\n",
    "        #                       width=800, height=475)\n",
    "        #return map_tiles * taxi_trips"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#TEST reproject dynamically with bounding box\n",
    "#!pip install paramnb\n",
    "#!conda install -y -c ioam parambokeh\n",
    "import param\n",
    "import paramnb\n",
    "#import parambokeh\n",
    "\n",
    "minAlt = round(df['z'].min().compute(), -2)\n",
    "maxAlt = round(df['z'].max().compute(), -2)\n",
    "minTrk = int(df['k'].min().compute())\n",
    "maxTrk = int(df['k'].max().compute())\n",
    "\n",
    "class IceSatExplorer(hv.streams.Stream):\n",
    "    #plotdims   = param.ObjectSelector(default=['x', 'y'], objects=[['x', 'y'], ['y', 'z'], ['x', 'z']])\n",
    "    colormap   = param.ObjectSelector(default=cm[\"fire\"], objects=cm.values())\n",
    "    altitude   = param.Range(default=(minAlt, maxAlt), bounds=(minAlt, maxAlt), doc=\"\"\"Elevation of ICESAT laser point\"\"\")\n",
    "    trackrange = param.Range(default=(minTrk, maxTrk), bounds=(minTrk, maxTrk), doc=\"\"\"ICESAT Track subset\"\"\" )\n",
    "    timerange  = param.Range()\n",
    "    \n",
    "    def make_view(self, plotdims=['x', 'y'], x_range=None, y_range=None, **kwargs):\n",
    "        datas = projected_gv.select(z=self.altitude, k=self.trackrange).to(hv.Dataset)\n",
    "        \n",
    "        hvmap = datas.to(hv.Points, kdims=plotdims, vdims=['i', 'k'], groupby=[])\n",
    "        #dsmap = datashade(hvmap, cmap=self.colormap, dynamic=False, x_range=x_range, y_range=y_range)\n",
    "        #gv_options = {'bgcolor':'black', 'show_grid':True}\n",
    "        #gvmap = dynspread(dsmap.opts(plot=gv_options))\n",
    "        gvmap = hvmap\n",
    "        return gvmap\n",
    "    \n",
    "    def z_y_views(self, x_range, y_range):\n",
    "        x_min = x_range[0]; x_max = x_range[1]\n",
    "        y_min = y_range[0]; y_max = y_range[1]\n",
    "        \n",
    "        datas = projected_gv.select(x=x_range, y=y_range, z=self.altitude, k=self.trackrange).to(hv.Dataset)\n",
    "        hvmap = datas.to(hv.Points, kdims=self.plotdims, vdims=['i', 'k'], groupby=[])\n",
    "        dsmap = datashade(hvmap, cmap=self.colormap, dynamic=False, x_range=x_range, y_range=y_range)\n",
    "        gv_options = {'bgcolor':'black', 'show_grid':True}\n",
    "        gvmap = dynspread(dsmap.opts(plot=gv_options))\n",
    "        \n",
    "        return gvmap\n",
    "        #return pts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%output size=150  #set output size, e.g. 200 = 2x the default output size\n",
    "#%%opts Points [tools=['box_select', 'lasso_select']]\n",
    "explorer = IceSatExplorer() \n",
    "paramnb.Widgets(explorer, callback=explorer.event)\n",
    "\n",
    "def z_y_points(x_range, y_range):\n",
    "    x_min = x_range[0]; x_max = x_range[1]\n",
    "    #y_min = y_range[0]; y_max = y_range[1]\n",
    "    pts = projected_gv.select(x=x_range).to(hv.Dataset)\n",
    "    pts = pts.to(hv.Points, kdims=['z', 'y'], groupby=[])\n",
    "    return pts\n",
    "\n",
    "x_y = explorer.make_view(plotdims=['x', 'y'])\n",
    "dmap = hv.DynamicMap(z_y_points, streams=[hv.streams.RangeXY(x_range=(-2500000, 2500000), source=x_y)], kdims=[])\n",
    "\n",
    "dynspread(datashade(x_y)) + dynspread(datashade(dmap))\n",
    "\n",
    "#xy_map = hv.DynamicMap(explorer.make_view, streams=[explorer, RangeXY()])\n",
    "#xy_map\n",
    "#hv.help(xy_map)\n",
    "\n",
    "#+ \\\n",
    "#datashade(dmap)\n",
    "#hv.DynamicMap(IceSatExplorer(plotdims=['x', 'z']).make_view, streams=[explorer, RangeXY(source=explorer.make_view)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raise ValueError('temp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matplotlib 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(pdData.loc[:,['x']], pdData.loc[:,['y']]); #2d plot, will show Greenland on top right and Antarctica at the bottom\n",
    "plt.show()\n",
    "plt.scatter(df.loc[:,['x']], df.loc[:,['y']]); #2d plot, showing data for Antarctica\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matplotlib 3D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.set_xlabel('X Label')\n",
    "ax.set_ylabel('Y Label')\n",
    "ax.set_zlabel('Z Label')\n",
    "mi=10000\n",
    "x=df.loc[:mi,['x']]\n",
    "y=df.loc[:mi,['y']]\n",
    "z=df.loc[:mi,['z']]\n",
    "ax.scatter(x, y, z, c=z, cmap=plt.cm.jet, zdir='z', s=2)\n",
    "fig;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%% Matplotlib 3D surface plot https://gis.stackexchange.com/a/66440/78212\n",
    "#ax = fig.add_subplot(111, projection='3d')\n",
    "#ax.plot_surface(x,y,z)\n",
    "#ax.contour(x,y,z)\n",
    "#fig;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Folium"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip3 install folium\n",
    "import folium\n",
    "print(folium.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#%%timeit\n",
    "map_osm = folium.Map(location=[-77.84651, 166.75710], zoom_start=3)\n",
    "#df[:1000].apply(lambda row:folium.CircleMarker(location=[row['y'], row[\"x\"]], popup=str(row['k']), radius=5).add_to(map_osm), axis=1)\n",
    "map_osm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1.52 s ± 6.52 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### laspy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip3 install laspy\n",
    "import laspy\n",
    "print(laspy.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "header = laspy.header.Header()\n",
    "outfile = laspy.file.File(os.environ['HOME']+\"/code/scripts/output.las\", mode=\"w\", header=header)\n",
    "\n",
    "allx = df.loc[:,['x']].values.flatten()\n",
    "ally = df.loc[:,['y']].values.flatten()\n",
    "allz = df.loc[:,['z']].values.flatten()\n",
    "\n",
    "xmin = np.floor(np.min(allx))\n",
    "ymin = np.floor(np.min(ally))\n",
    "zmin = np.floor(np.min(allz))\n",
    "\n",
    "outfile.header.offset = [xmin,ymin,zmin]\n",
    "outfile.header.scale = [0.001,0.001,0.001]\n",
    "\n",
    "outfile.x = allx\n",
    "outfile.y = ally\n",
    "outfile.z = allz\n",
    "\n",
    "outfile.close()"
   ]
  },
  {
   "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": [
    "# %%\n",
    "latvar = f['/Data_1HZ/Geolocation/d_lat']\n",
    "latitude = latvar[:]\n",
    "lat_vr = [latvar.attrs['valid_min'], latvar.attrs['valid_max']]\n",
    "\n",
    "lonvar = f['/Data_1HZ/Geolocation/d_lon']\n",
    "longitude = lonvar[:]\n",
    "lon_vr = [lonvar.attrs['valid_min'], lonvar.attrs['valid_max']]\n",
    "\n",
    "latlon = np.dstack((latvar, lonvar))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def h5_treeview(h5file):\n",
    "        assert isinstance(h5file, (h5py._hl.files.File))\n",
    "\n",
    "        for k1 in sorted(h5file.keys()):\n",
    "            print(k1)\n",
    "            for k2 in f[k1]:\n",
    "                print(\"'--\", k2)\n",
    "\n",
    "h5_treeview(f)\n",
    "f['ANCILLARY_DATA']\n",
    "\n",
    "f[datagroup].visititems(lambda name, object: print(name.count('/')*\"    \"+\"'--\"+name, object))\n",
    "f[datagroup].visititems(lambda name, object: (None))"
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}