Mars_5Bins_EqualArea_by_Elevations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate 5 equal area bins of Mars elevation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using equal area pixels (from a Sinusoidal projection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "read from MOLA DEM into numpy array "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ds = gdal.Open(\"megt90n000_sinu_topo.tif\")\n",
    "elev = np.array(ds.GetRasterBand(1).ReadAsArray())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#tested reading table also\n",
    "#csv = np.genfromtxt ('mars_16ppd_sinu_table_ascii.csv', comments=\"p\", delimiter=\",\")\n",
    "#elev = csv[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "#remove NoData from file\n",
    "elev = elev[elev != 32767]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "elev = elev.flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10560695"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "elev.size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-8177, -8155, -8149, ..., 21162, 21164, 21171], dtype=int16)"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "elev_sort = np.sort(elev)\n",
    "elev_sort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#function simply breaks into like \"equal areas\" sizes because of map projection\n",
    "def equal_bin(N, m):\n",
    "    sep = (N.size/float(m))*np.arange(1,m+1)\n",
    "    idx = sep.searchsorted(np.arange(N.size))\n",
    "    return idx[N.argsort().argsort()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2112140, 2112139, 2112139, 2112139, 2112138], dtype=int64)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# setup the index into the 5 bins\n",
    "bins = equal_bin(elev, 5)\n",
    "bin_index = np.bincount(bins)\n",
    "bin_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2112139, 4224278, 6336417, 8448556, 10560694)"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Just add up index maxes as a look-up for the sorted elevation array\n",
    "bin1_imax = bin_index[0] - 1\n",
    "bin2_imax = bin_index[1] + bin1_imax\n",
    "bin3_imax = bin_index[2] + bin2_imax\n",
    "bin4_imax = bin_index[3] + bin3_imax\n",
    "bin5_imax = bin_index[4] + bin4_imax\n",
    "\n",
    "bin1_imax, bin2_imax, bin3_imax, bin4_imax, bin5_imax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "bin1_min, bin1_max = elev_sort[0],           elev_sort[bin1_imax]\n",
    "bin2_min, bin2_max = elev_sort[bin1_imax+1], elev_sort[bin2_imax]\n",
    "bin3_min, bin3_max = elev_sort[bin2_imax+1], elev_sort[bin3_imax]\n",
    "bin4_min, bin4_max = elev_sort[bin3_imax+1], elev_sort[bin4_imax]\n",
    "bin5_min, bin5_max = elev_sort[bin4_imax+1], elev_sort[bin5_imax]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-8177 -3702\n",
      "-3701 -1475\n",
      "-1474 794\n",
      "795 1819\n",
      "1820 21171\n"
     ]
    }
   ],
   "source": [
    "print(bin1_min, bin1_max)\n",
    "print(bin2_min+1, bin2_max)\n",
    "print(bin3_min+1, bin3_max)\n",
    "print(bin4_min+1, bin4_max)\n",
    "print(bin5_min+1, bin5_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Done, now apply elevation ranges to reclassify image and colorize in a GIS"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}