{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# 2D histogram of the optical depth $\\tau$\n", "\n", "Below I calculate the 2-d and averaged 1-d spectra for the optical depth, which gives the penetration\n", "depth of photons through a cloud, and is closely related to cloud thickness" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\",category=FutureWarning)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "download a17.nc: size is 17 Mbytes\n" ] } ], "source": [ "from matplotlib import pyplot as plt\n", "import urllib\n", "import os\n", "filelist=['a17.nc']\n", "data_download=True\n", "if data_download:\n", " for the_file in filelist:\n", " url='http://clouds.eos.ubc.ca/~phil/docs/atsc500/data/{}'.format(the_file)\n", " urllib.request.urlretrieve(url,the_file)\n", "print(\"download {}: size is {:6.2g} Mbytes\".format(the_file,os.path.getsize(the_file)*1.e-6))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'netCDF4'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnetCDF4\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mDataset\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mDataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilelist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnc\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0mtau\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'tau'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m...\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'netCDF4'" ] } ], "source": [ "from netCDF4 import Dataset\n", "with Dataset(filelist[0]) as nc:\n", " tau=nc.variables['tau'][...]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Character of the optical depth field\n", "\n", "The image below shows one of the marine boundary layer landsat scenes analyzed in \n", "[Lewis et al., 2004](http://onlinelibrary.wiley.com/doi/10.1029/2003JD003742/full)\n", "\n", "It is a 2048 x 2048 pixel image taken by Landsat 7, with the visible reflectivity converted to\n", "cloud optical depth. The pixels are 25 m x 25 m, so the scene extends for about 50 km x 50 km" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'tau' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m13\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'landsat a17'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mim0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7\u001b[0m \u001b[0mim1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtau\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'histogram of tau values'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'tau' is not defined" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwIAAAGrCAYAAABzMCC3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAXxUlEQVR4nO3df4zkd33f8debM1b4lUDjA4FtGicyGKvCKRw/GpXglKbYJK0VBVF+FBIU1XIT0qhqWqyoDa1QG1K1FUIYHINcSlvFjYKVGOJgpUqBSI5bnyswGAq9mGBf7JQzv39EoDPv/jFjuixr79x5d2/n3o+HdNLOfL+789mP7uZ9z53ZmeruAAAAszziVC8AAADYe0IAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICHAKVVVf1pVf3MXvu4PVFVX1Rk7/bUBAE4HQgAewsMNlap6WVXdXFVfr6oPbDr2gqr66qY/XVU//bAXDgCwDSEAu+vzSd6c5E2bD3T3H3X3Yx/4k+Qnk3w1yfv3eI0AwEBCgH2jqp5bVX9cVV+sqnur6q1VdeaG411VV1TV/6mqL1TVVVVVy2MHqurfVtV9VXVnkp/Y9LV/tqrurKqvVNWnq+pVy+t/qKr+sKo+t/zc/1JVj18e+09Jnprkvcuf1v/TLdb8hKp6X1UdW67pfVV1zgPHu/u/dfdvJblnhS34mSS/3d1fO/HdAwA4MUKA/eT+JP8oyVlJ/lqSFyX5+U3n/GSS5yS5KMnLkrx4ef3fXx77q0kOJXnpA59QVY9J8pYkl3b345L8SJIPP3A4ya8leUqSZyQ5N8m/SJLufnWSu5L87eVP7f/NFmt+RJL/kOQvZxENf5HkrSf6jVfVo5dr/o8n+rkAACdDCLBvdPdt3X1Ldx/v7j9N8htJXrjptDd19xe7+64k/z3JDy+vf1mSN3f33d39+Sz+c7/Rt5L8lap6VHff2913LG/zSHf/QXd/o7uPJfn3W9zmQ635c939nu7+end/Jcm/OpHP3+Cnk9yX5IMn8bkAACdMCLBvVNXTlk+t+fOq+nKSf53FowMb/fmGj7+e5LHLj5+S5O4Nxz7zwAfLp9r83SRXJLm3qn6vqi5Y3uYTq+q6qvqz5W3+5y1u86HW/Oiq+o2q+szy8z+U5PFVdWDVr7H0M0ne3d19gp8HAHBShAD7yduT/O8k53f39yb5lSyeurOKe7N4Ws8DnrrxYHff1N0/nuTJy9t4x/LQryXpJM9c3ubf23Sb2/3H/B8neXqS5y0//0eX16+67lTVuUkuTvLuVT8HAODhEgLsJ49L8uUkX13+xP4fnMDn/laSf1hV51TVE5Jc+cCBqnpSVf2d5e8KfCOLV+a5f8NtfjXJF6vq7CT/ZNPX/b9JfnCbNf/F8vP/UpI3bDy4/CXm70lyRpJHVNX3VNUjN32NVye5ubv/5AS+XwCAh0UIsJ/8cpJXJvlKFj+x/68n8LnvSHJTko8k+V9Jrt9w7BFZ/OT+nixezvOF+f+/hPwvkzwryZeS/N6mz0sWjxj8s+UrGf3yFrf75iSPyuL5/bfku1/689VZhMLbk7xg+fE7Np3zmvglYQBgj5WnJAMAwDweEQAAgIG2DYGquraqPltVH3uQ41VVb6mqI1V1e1U9a+eXCcB+Z14ArJdVHhF4V5JLHuL4pUnOX/65PIvnQgMwz7tiXgCsjW1DoLs/lMUvWD6Yy7J8/fPuviWL11B/8k4tEID1YF4ArJczduBrnJ3vfCOno8vr7t18YlVdnsVPgfKYxzzm2RdccMEO3DzA6em22267r7sPnup17CDzAmAXnOy82IkQ2OqNk7Z8KaLuvibJNUly6NChPnz48A7cPMDpqao+s/1Za8W8ANgFJzsvduJVg47mO9/R9ZwsXq8dADYyLwD2kZ0IgRuSvGb5ahDPT/Kl7v6uh3kBGM+8ANhHtn1qUFX9ZpKLk5xVVUeTvCHJI5Oku69OcmOSlyQ5kuTrSV67W4sFYP8yLwDWy7Yh0N2v2OZ4J/mFHVsRAGvJvABYL95ZGAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIFWCoGquqSqPllVR6rqyi2Of19VvbeqPlJVd1TVa3d+qQDsZ2YFwHrZNgSq6kCSq5JcmuTCJK+oqgs3nfYLST7e3RcluTjJv6uqM3d4rQDsU2YFwPpZ5RGB5yY50t13dvc3k1yX5LJN53SSx1VVJXlsks8nOb6jKwVgPzMrANbMKiFwdpK7N1w+urxuo7cmeUaSe5J8NMkvdfe3Nn+hqrq8qg5X1eFjx46d5JIB2Id2bFYk5gXAXlglBGqL63rT5Rcn+XCSpyT54SRvrarv/a5P6r6muw9196GDBw+e8GIB2Ld2bFYk5gXAXlglBI4mOXfD5XOy+GnORq9Ncn0vHEny6SQX7MwSAVgDZgXAmlklBG5Ncn5Vnbf8pa6XJ7lh0zl3JXlRklTVk5I8PcmdO7lQAPY1swJgzZyx3QndfbyqXpfkpiQHklzb3XdU1RXL41cneWOSd1XVR7N4ePj13X3fLq4bgH3ErABYP9uGQJJ0941Jbtx03dUbPr4nyd/a2aUBsE7MCoD14p2FAQBgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGGilEKiqS6rqk1V1pKqufJBzLq6qD1fVHVX1wZ1dJgD7nVkBsF7O2O6EqjqQ5KokP57kaJJbq+qG7v74hnMen+RtSS7p7ruq6om7tWAA9h+zAmD9rPKIwHOTHOnuO7v7m0muS3LZpnNemeT67r4rSbr7szu7TAD2ObMCYM2sEgJnJ7l7w+Wjy+s2elqSJ1TVB6rqtqp6zVZfqKour6rDVXX42LFjJ7diAPajHZsViXkBsBdWCYHa4rredPmMJM9O8hNJXpzkn1fV077rk7qv6e5D3X3o4MGDJ7xYAPatHZsViXkBsBe2/R2BLH6qc+6Gy+ckuWeLc+7r7q8l+VpVfSjJRUk+tSOrBGC/MysA1swqjwjcmuT8qjqvqs5M8vIkN2w653eTvKCqzqiqRyd5XpJP7OxSAdjHzAqANbPtIwLdfbyqXpfkpiQHklzb3XdU1RXL41d39yeq6v1Jbk/yrSTv7O6P7ebCAdg/zAqA9VPdm5/CuTcOHTrUhw8fPiW3DbAOquq27j50qtdxqpkXAA/tZOeFdxYGAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADrRQCVXVJVX2yqo5U1ZUPcd5zqur+qnrpzi0RgHVgVgCsl21DoKoOJLkqyaVJLkzyiqq68EHO+/UkN+30IgHY38wKgPWzyiMCz01ypLvv7O5vJrkuyWVbnPeLSd6T5LM7uD4A1oNZAbBmVgmBs5PcveHy0eV131ZVZyf5qSRXP9QXqqrLq+pwVR0+duzYia4VgP1rx2bF8lzzAmCXrRICtcV1venym5O8vrvvf6gv1N3XdPeh7j508ODBVdcIwP63Y7MiMS8A9sIZK5xzNMm5Gy6fk+SeTeccSnJdVSXJWUleUlXHu/t3dmSVAOx3ZgXAmlklBG5Ncn5VnZfkz5K8PMkrN57Q3ec98HFVvSvJ+9yxA4xiVgCsmW1DoLuPV9XrsniFhwNJru3uO6rqiuXxbZ/rCcDpzawAWD+rPCKQ7r4xyY2brtvyTr27f/bhLwuAdWNWAKwX7ywMAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQCuFQFVdUlWfrKojVXXlFsdfVVW3L//cXFUX7fxSAdjPzAqA9bJtCFTVgSRXJbk0yYVJXlFVF2467dNJXtjdz0zyxiTX7PRCAdi/zAqA9bPKIwLPTXKku+/s7m8muS7JZRtP6O6bu/sLy4u3JDlnZ5cJwD5nVgCsmVVC4Owkd2+4fHR53YP5uSS/v9WBqrq8qg5X1eFjx46tvkoA9rsdmxWJeQGwF1YJgdriut7yxKofy+LO/fVbHe/ua7r7UHcfOnjw4OqrBGC/27FZkZgXAHvhjBXOOZrk3A2Xz0lyz+aTquqZSd6Z5NLu/tzOLA+ANWFWAKyZVR4RuDXJ+VV1XlWdmeTlSW7YeEJVPTXJ9Ule3d2f2vllArDPmRUAa2bbRwS6+3hVvS7JTUkOJLm2u++oqiuWx69O8qtJvj/J26oqSY5396HdWzYA+4lZAbB+qnvLp3DuukOHDvXhw4dPyW0DrIOqus1/lM0LgO2c7LzwzsIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMJAQAAGAgIQAAAAMJAQAAGEgIAADAQEIAAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABhICAAAwEBCAAAABhICAAAwkBAAAICBhAAAAAwkBAAAYCAhAAAAAwkBAAAYSAgAAMBAQgAAAAYSAgAAMJAQAACAgYQAAAAMtFIIVNUlVfXJqjpSVVducbyq6i3L47dX1bN2fqkA7GdmBcB62TYEqupAkquSXJrkwiSvqKoLN512aZLzl38uT/L2HV4nAPuYWQGwflZ5ROC5SY50953d/c0k1yW5bNM5lyV5dy/ckuTxVfXkHV4rAPuXWQGwZs5Y4Zyzk9y94fLRJM9b4Zyzk9y78aSqujyLnwIlyTeq6mMntNrT01lJ7jvVi9gH7MOCfViwDwtPP9ULOAE7NisS8+JB+HexYB8W7MOCfVg4qXmxSgjUFtf1SZyT7r4myTVJUlWHu/vQCrd/WrMPC/ZhwT4s2IeFqjp8qtdwAnZsViTmxVbsw4J9WLAPC/Zh4WTnxSpPDTqa5NwNl89Jcs9JnAPA6cusAFgzq4TArUnOr6rzqurMJC9PcsOmc25I8prlK0I8P8mXuvu7HuoF4LRlVgCsmW2fGtTdx6vqdUluSnIgybXdfUdVXbE8fnWSG5O8JMmRJF9P8toVbvuak1716cU+LNiHBfuwYB8W1mYfdnFWJGu0D7vMPizYhwX7sGAfFk5qH6p7y6dnAgAApzHvLAwAAAMJAQAAGGjXQ8Bbzi+ssA+vWn7/t1fVzVV10alY527bbh82nPecqrq/ql66l+vbK6vsQ1VdXFUfrqo7quqDe73GvbDCv4vvq6r3VtVHlvuw6nPK10ZVXVtVn32w18l3H/nt4/YhZsUW55kVZsWIWZHs0rzo7l37k8UvjP1Jkh9McmaSjyS5cNM5L0ny+1m8vvTzk/yP3VzTqfiz4j78SJInLD++dOo+bDjvD7P4xcKXnup1n6K/D49P8vEkT11efuKpXvcp2odfSfLry48PJvl8kjNP9dp3eB9+NMmzknzsQY67j7QPG88xK77zPLPCrBgxK5bf247Pi91+RMBbzi9suw/dfXN3f2F58ZYsXl/7dLPK34ck+cUk70ny2b1c3B5aZR9emeT67r4rSbr7dNyLVfahkzyuqirJY7O4cz++t8vcXd39oSy+rwfjPnLBPsSs2MSsMCsecNrPimR35sVuh8CDvZ38iZ6z7k70e/y5LIrudLPtPlTV2Ul+KsnVe7iuvbbK34enJXlCVX2gqm6rqtfs2er2zir78NYkz8jiTac+muSXuvtbe7O8fcN95OrnrDuzYsGsWDArFsyK1Z3w/eS27yPwMO3oW86vsZW/x6r6sSzu3P/6rq7o1FhlH96c5PXdff8i7E9Lq+zDGUmeneRFSR6V5I+r6pbu/tRuL24PrbIPL07y4SR/I8kPJfmDqvqj7v7ybi9uH3Efufo5686sWDArFsyKBbNidSd8P7nbIeAt5xdW+h6r6plJ3pnk0u7+3B6tbS+tsg+Hkly3vGM/K8lLqup4d//O3ixxT6z67+K+7v5akq9V1YeSXJTkdLpzX2UfXpvkTb148uORqvp0kguS/M+9WeK+4D5y9XPWnVmxYFYsmBULZsXqTvh+crefGuQt5xe23YeqemqS65O8+jQr+Y223YfuPq+7f6C7fyDJbyf5+dPsjj1Z7d/F7yZ5QVWdUVWPTvK8JJ/Y43XutlX24a4sftKVqnpSkqcnuXNPV3nquY9csA8xKx5gVnybWbFgViyc8P3krj4i0Lv7lvNrY8V9+NUk35/kbcufcBzv7kOnas27YcV9OO2tsg/d/Ymqen+S25N8K8k7u3vLlwtbVyv+fXhjkndV1UezeMjz9d193ylb9C6oqt9McnGSs6rqaJI3JHlk4j7SrDArYlaYFWbFt+3GvKjFoygAAMAk3lkYAAAGEgIAADCQEAAAgIGEAAAADCQEAABgICEAAAADCQEAABjo/wH4/wMOtEZe9AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", "plt.close('all')\n", "fig,ax=plt.subplots(1,2,figsize=(13,7))\n", "ax[0].set_title('landsat a17')\n", "im0=ax[0].imshow(tau)\n", "im1=ax[1].hist(tau.ravel())\n", "ax[1].set_title('histogram of tau values')\n", "divider = make_axes_locatable(ax[0])\n", "cax = divider.append_axes(\"bottom\", size=\"5%\", pad=0.35)\n", "out=fig.colorbar(im0,orientation='horizontal',cax=cax)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## ubc_fft class" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the next cell I define a class that calculates the 2-d fft for a square image\n", "\n", "in the method ```power_spectrum``` we calculate both the 2d fft and the power spectrum\n", "and save them as class attributes. In the method ```annular_average``` I take the power spectrum,\n", "which is the two-dimensional field $E(k_x, k_y)$ (in cartesian coordinates) or $E(k,\\theta)$ (in polar coordinates).\n", "In the method ```annular_avg``` I take the average\n", "\n", "$$\n", "\\overline{E}(k) = \\int_0^{2\\pi} E(k, \\theta) d\\theta\n", "$$\n", "\n", "and plot that average with the method ```graph_spectrum```" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true, "lines_to_next_cell": 2 }, "outputs": [ { "ename": "ModuleNotFoundError", "evalue": "No module named 'netCDF4'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mModuleNotFoundError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mnetCDF4\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mDataset\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mmath\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mfft\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mmatplotlib\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpyplot\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'netCDF4'" ] } ], "source": [ "from netCDF4 import Dataset\n", "import numpy as np\n", "import math\n", "from numpy import fft\n", "from matplotlib import pyplot as plt\n", "\n", "\n", "class ubc_fft:\n", "\n", " def __init__(self, filename, var, scale):\n", " \"\"\"\n", " Input filename, var=variable name, \n", " scale= the size of the pixel in km\n", "\n", " Constructer opens the netcdf file, reads the data and\n", " saves the twodimensional fft\n", " \"\"\"\n", " with Dataset(filename,'r') as fin:\n", " data = fin.variables[var][...]\n", " data = data - data.mean()\n", " if data.shape[0] != data.shape[1]:\n", " raise ValueError('expecting square matrix')\n", " self.xdim = data.shape[0] # size of each row of the array\n", " self.midpoint = int(math.floor(self.xdim/2))\n", " root,suffix = filename.split('.')\n", " self.filename = root\n", " self.var = var\n", " self.scale = float(scale)\n", " self.data = data\n", " self.fft_data = fft.fft2(self.data)\n", " \n", " def power_spectrum(self):\n", " \"\"\"\n", " calculate the power spectrum for the 2-dimensional field\n", " \"\"\"\n", " #\n", " # fft_shift moves the zero frequency point to the middle\n", " # of the array \n", " #\n", " fft_shift = fft.fftshift(self.fft_data)\n", " spectral_dens = fft_shift*np.conjugate(fft_shift)/(self.xdim*self.xdim)\n", " spectral_dens = spectral_dens.real\n", " #\n", " # dimensional wavenumbers for 2dim spectrum (need only the kx\n", " # dimensional since image is square\n", " #\n", " k_vals = np.arange(0,(self.midpoint))+1\n", " k_vals = (k_vals-self.midpoint)/(self.xdim*self.scale)\n", " self.spectral_dens=spectral_dens\n", " self.k_vals=k_vals\n", "\n", " def annular_avg(self,avg_binwidth):\n", " \"\"\" \n", " integrate the 2-d power spectrum around a series of rings \n", " of radius kradial and average into a set of 1-dimensional\n", " radial bins\n", " \"\"\"\n", " #\n", " # define the k axis which is the radius in the 2-d polar version of E\n", " #\n", " numbins = int(round((math.sqrt(2)*self.xdim/avg_binwidth),0)+1)\n", "\n", " avg_spec = np.zeros(numbins,np.float64)\n", " bin_count = np.zeros(numbins,np.float64)\n", "\n", " print(\"\\t- INTEGRATING... \")\n", " for i in range(self.xdim):\n", " if (i%100) == 0:\n", " print(\"\\t\\trow: {} completed\".format(i))\n", " for j in range(self.xdim):\n", " kradial = math.sqrt(((i+1)-self.xdim/2)**2+((j+1)-self.xdim/2)**2)\n", " bin_num = int(math.floor(kradial/avg_binwidth))\n", " avg_spec[bin_num]=avg_spec[bin_num]+ kradial*self.spectral_dens[i,j]\n", " bin_count[bin_num]+=1\n", "\n", " for i in range(numbins):\n", " if bin_count[i]>0:\n", " avg_spec[i]=avg_spec[i]*avg_binwidth/bin_count[i]/(4*(math.pi**2))\n", " self.avg_spec=avg_spec\n", " #\n", " # dimensional wavenumbers for 1-d average spectrum\n", " #\n", " self.k_bins=np.arange(numbins)+1\n", " self.k_bins = self.k_bins[0:self.midpoint]\n", " self.avg_spec = self.avg_spec[0:self.midpoint]\n", "\n", " \n", " \n", " def graph_spectrum(self, kol_slope=-5./3., kol_offset=1., \\\n", " title=None):\n", " \"\"\"\n", " graph the annular average and compare it to Kolmogorov -5/3\n", " \"\"\"\n", " avg_spec=self.avg_spec\n", " delta_k = 1./self.scale # 1./km (1/0.025 for landsat 25 meter pixels)\n", " nyquist = delta_k * 0.5\n", " knum = self.k_bins * (nyquist/float(len(self.k_bins)))# k = w/(25m)\n", " #\n", " # draw the -5/3 line through a give spot\n", " #\n", " kol = kol_offset*(knum**kol_slope)\n", " fig,ax=plt.subplots(1,1,figsize=(8,8))\n", " ax.loglog(knum,avg_spec,'r-',label='power')\n", " ax.loglog(knum,kol,'k-',label=\"$k^{-5/3}$\")\n", " ax.set(title=title,xlabel='k (1/km)',ylabel='$E_k$')\n", " ax.legend()\n", " self.plotax=ax" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'ubc_fft' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'all'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstyle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0muse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ggplot'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0moutput\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mubc_fft\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'a17.nc'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'tau'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0.025\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpower_spectrum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'ubc_fft' is not defined" ] } ], "source": [ "plt.close('all')\n", "plt.style.use('ggplot')\n", "output = ubc_fft('a17.nc','tau',0.025)\n", "output.power_spectrum()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'np' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m7\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'landsat a17'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mim0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mspectral_dens\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_title\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'log10 of the 2-d power spectrum'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mdivider\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmake_axes_locatable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbQAAAGsCAYAAAC1sVKwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAaU0lEQVR4nO3df2xV9f3H8ddtb6CBllruTVsKCFKQH9kMhEpJJ4xKU5gGbUaCbGQZY41oZYKTRUDAimNrHMivgMPQlOnc9gduw2STkBth/KiTsrYEMAzKEMZawN6CUH7ae8/3j682Xlu4lwP9wbvPR0Li6f303A9vhaf33MvB4ziOIwAA7nFxHb0BAADuBoIGADCBoAEATCBoAAATCBoAwASCBgAwgaDBtJkzZyovL69dnqu4uFiDBw9ul+cC0BJBAzqpPXv2yOPx6NNPP70r55s3b56ys7PVo0cPeb3eFo/v3LlTHo+n1R+/+c1v7soegLZE0IAuIhQK6Yc//KGKiopafTwnJ0d1dXURP9544w3FxcVp2rRp7bxb4PYRNHQplZWV+t73vqfU1FQlJibq4Ycf1rZt2yLWDBw4UEuXLtXcuXPVu3dvpaWlaf78+QqFQs1rrl+/rmeffVbJyclKSUnRs88+q+vXr0ec5/Dhw5o0aZLuu+8+9ezZU8OHD9c777zT/PiaNWs0cuRIJSYmKj09XdOnT1ddXZ0k6dNPP9W4ceMkSQ888IA8Ho8mTJhw05/Xrc71lXXr1mnu3Ln61re+1eo5unXrpvT09IgfW7Zs0eTJkzVgwIDowwU6GEFDl3Lx4kVNnz5dO3fuVGVlpSZNmqQnnnhCR48ejVi3bt069enTRx9//LHWrl2r1atX6+23325+fMGCBXrvvff09ttv66OPPlLPnj21fv36iHP84Ac/kM/nU3l5uQ4ePKg33nhDKSkpEWtWrFihgwcP6i9/+YtOnTql6dOnS5L69++vrVu3SpL27dunuro6/fnPf77lz+1m53Lr4MGDKi8v1+zZs+/oPEC7cQDDfvzjHzsTJ0685ZqHHnrI+eUvf9l8PGDAAGfKlCkRayZNmuRMnz7dcRzHaWxsdLp37+689dZbEWtGjx7tZGZmNh/36tXLKSsri3mvlZWVjiTn9OnTjuM4zu7dux1JzokTJ2I+x83O9XVlZWVOfHx81HM899xzTt++fZ2mpqbbfn6gI/AKDV3KZ599pqKiIg0bNkz33XefEhMTdfjwYZ08eTJi3ciRIyOO+/btq7Nnz0qSjh8/ruvXrysnJydizSOPPBJxPH/+fBUWFmrChAkqLi5WZWVlxOM7d+7UpEmT1L9/fyUlJTV//zf3Eou7eS5JunLlin7/+9+rsLBQ8fHxrs4BtDeChi5l5syZ2r17t15//XXt3r1b1dXVGjlypG7cuBGxrlu3bhHHHo9H4XBYkuR8+RdUeDyeWz7XkiVLdPToUU2bNk2HDh3S2LFjtXjxYknSqVOn9Nhjj2ngwIH605/+pP379+v999+XpBZ7ieZunusrf/zjH9XY2KjCwkJX3w90BIKGLmXXrl0qKirSE088oW9/+9vq06eP/vOf/9zWOQYPHqxu3bpp7969EV8vLy9vsXbQoEEqKirSli1btGzZMr355puSpIqKCl29elWrV6/Wd77zHQ0dOrT5FeBXvorq1z+M0ppYznW7Nm7cqMcff1z9+vW7o/MA7anlH0YBDBs6dKjeffddPfLIIwqFQlq6dGnUYHxTz5499cwzz2jx4sVKS0vT0KFDVVpaqiNHjig1NVWS1NjYqJdeeklTp07VAw88oAsXLmjbtm0aMWKEJGnIkCHyeDxauXKlZsyYoQMHDmjZsmURzzNgwADFxcXp73//u5566il1795dycnJLfYTy7kkqaamRo2NjTp16pQkqbq6WtL/BzoxMbF5XVVVlSoqKvS3v/3ttuYCdDReoaFLKSsrUzgc1pgxY1RQUKDJkyfr4Ycfvu3zlJSUqKCgQD/60Y80ZswYXbhwQc8991zz416vV+fPn9dPf/pTDR8+XJMmTVJaWpr+8Ic/SJIeeughrVu3Ths3btSIESO0YsUKrV69OuI50tLS9Otf/1olJSXq06ePnnzyyVb3Esu5JKmwsFCjRo3SK6+8olAopFGjRmnUqFHav39/xLqNGzfq/vvv1+TJk297LkBH8jgOf2M1AODexys0AIAJUd9D27BhgyorK5WcnKyVK1e2eNxxHJWVlamqqkrdu3dXUVGRBg0a1CabBQDgZqK+QpswYYIWLVp008erqqp05swZrV27Vk8//bQ2bdp0VzcIAEAsogZtxIgREZ+A+qb9+/dr/Pjx8ng8evDBB3X58mWdP3/+rm4SAIBo7vhj+w0NDfL7/c3HPp9PDQ0NLe5ZJ0mBQECBQEDS/39KDACAu+WOg9bahyRvdgeFvLy8iL9ssba29k6fvsvx+/2qr6/v6G3cc5ibe8zOHebmTkZGhuvvveNPOfp8voh/acFgsNVXZwAAtKU7DlpWVpZ27dolx3F09OhR9ejRg6ABANpd1EuOq1ev1ieffKJLly7pmWee0bRp09TU1CRJys/P16hRo1RZWannn39e3bp1u+nfhgsAQFuKGrR58+bd8nGPx8MduQEAHY47hQAATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEbyyLqqurVVZWpnA4rIkTJ6qgoCDi8StXrmjt2rUKBoMKhUKaMmWKcnNz22TDAAC0JmrQwuGwSktLtXjxYvl8Pi1cuFBZWVnq169f85pt27apX79+WrBggS5evKi5c+dq3Lhx8npj6iUAAHcs6iXHmpoapaenKy0tTV6vVzk5OaqoqIhY4/F4dO3aNTmOo2vXrikxMVFxcVzNBAC0n6gvoRoaGuTz+ZqPfT6fjh07FrFm8uTJev311zV79mxdvXpVL7zwQqtBCwQCCgQCkqSSkhL5/f473X+X4/V6mZsLzM09ZucOc2t/UYPmOE6Lr3k8nojjAwcOaMCAAVq6dKnOnj2r1157TcOGDVOPHj0i1uXl5SkvL6/5uL6+3u2+uyy/38/cXGBu7jE7d5ibOxkZGa6/N+p1QZ/Pp2Aw2HwcDAaVkpISsWbHjh3Kzs6Wx+NRenq6UlNTVVtb63pTAADcrqhBy8zMVF1dnc6dO6empiaVl5crKysrYo3f79fBgwclSRcuXFBtba1SU1PbZscAALQi6iXH+Ph4zZo1S8uXL1c4HFZubq769++v7du3S5Ly8/M1depUbdiwQS+++KIkacaMGerVq1fb7hwAgK/xOK29SdZOuCx5+7gu7w5zc4/ZucPc3GnT99AAALgXEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYII3lkXV1dUqKytTOBzWxIkTVVBQ0GLN4cOHtXnzZoVCISUlJenVV1+965sFAOBmogYtHA6rtLRUixcvls/n08KFC5WVlaV+/fo1r7l8+bI2bdqkl19+WX6/X59//nmbbhoAgG+KesmxpqZG6enpSktLk9frVU5OjioqKiLW7NmzR9nZ2fL7/ZKk5OTkttktAAA3EfUVWkNDg3w+X/Oxz+fTsWPHItbU1dWpqalJxcXFunr1qh577DF997vfbXGuQCCgQCAgSSopKWkOIGLn9XqZmwvMzT1m5w5za39Rg+Y4TouveTyeiONQKKQTJ05oyZIlunHjhhYvXqwhQ4YoIyMjYl1eXp7y8vKaj+vr693uu8vy+/3MzQXm5h6zc4e5ufPNbtyOqEHz+XwKBoPNx8FgUCkpKS3WJCUlKSEhQQkJCRo+fLhOnjx5RxsDAOB2RH0PLTMzU3V1dTp37pyamppUXl6urKysiDVZWVk6cuSIQqGQrl+/rpqaGvXt27fNNg0AwDdFfYUWHx+vWbNmafny5QqHw8rNzVX//v21fft2SVJ+fr769eunkSNHav78+YqLi9Ojjz6q+++/v803DwDAVzxOa2+StZPa2tqOeup7Ftfl3WFu7jE7d5ibO3fyVhV3CgEAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYEJMQauurtbcuXP1s5/9TH/9619vuq6mpkZPPfWU/vnPf961DQIAEIuoQQuHwyotLdWiRYu0atUq7d27V6dPn2513bvvvquRI0e2yUYBALiVqEGrqalRenq60tLS5PV6lZOTo4qKihbrPvjgA2VnZ6tXr15tslEAAG7FG21BQ0ODfD5f87HP59OxY8darNm3b59eeeUVvfnmmzc9VyAQUCAQkCSVlJTI7/e73XeX5fV6mZsLzM09ZucOc2t/UYPmOE6Lr3k8nojjzZs3a8aMGYqLu/ULvry8POXl5TUf19fXx7pPfMnv9zM3F5ibe8zOHebmTkZGhuvvjRo0n8+nYDDYfBwMBpWSkhKx5vjx41qzZo0k6eLFi6qqqlJcXJzGjBnjemMAANyOqEHLzMxUXV2dzp07p969e6u8vFzPP/98xJr169dH/PPo0aOJGQCgXUUNWnx8vGbNmqXly5crHA4rNzdX/fv31/bt2yVJ+fn5bb5JAACi8TitvUnWTmprazvqqe9ZXJd3h7m5x+zcYW7u3Ml7aNwpBABgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACYQNACACQQNAGACQQMAmEDQAAAmEDQAgAkEDQBgAkEDAJhA0AAAJhA0AIAJBA0AYAJBAwCYQNAAACZ4Y1lUXV2tsrIyhcNhTZw4UQUFBRGP7969W1u3bpUkJSQkqLCwUAMHDrzrmwUA4GaivkILh8MqLS3VokWLtGrVKu3du1enT5+OWJOamqri4mKtWLFCU6dO1VtvvdVmGwYAoDVRg1ZTU6P09HSlpaXJ6/UqJydHFRUVEWuGDh2qxMRESdKQIUMUDAbbZrcAANxE1EuODQ0N8vl8zcc+n0/Hjh276foPP/xQo0aNavWxQCCgQCAgSSopKZHf77/d/XZ5Xq+XubnA3Nxjdu4wt/YXNWiO47T4msfjaXXtoUOHtGPHDi1btqzVx/Py8pSXl9d8XF9fH+s+8SW/38/cXGBu7jE7d5ibOxkZGa6/N+olR5/PF3EJMRgMKiUlpcW6kydPauPGjfrFL36hpKQk1xsCAMCNqEHLzMxUXV2dzp07p6amJpWXlysrKytiTX19vVasWKE5c+bcUV0BAHAr6iXH+Ph4zZo1S8uXL1c4HFZubq769++v7du3S5Ly8/O1ZcsWNTY2atOmTc3fU1JS0rY7BwDgazxOa2+StZPa2tqOeup7Ftfl3WFu7jE7d5ibO236HhoAAPcCggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAATCBoAwASCBgAwgaABAEwgaAAAEwgaAMAEggYAMIGgAQBMIGgAABMIGgDABIIGADCBoAEATPDGsqi6ulplZWUKh8OaOHGiCgoKIh53HEdlZWWqqqpS9+7dVVRUpEGDBrXJhgEAaE3UV2jhcFilpaVatGiRVq1apb179+r06dMRa6qqqnTmzBmtXbtWTz/9tDZt2tRmGwYAoDVRg1ZTU6P09HSlpaXJ6/UqJydHFRUVEWv279+v8ePHy+Px6MEHH9Tly5d1/vz5Nts0AADfFPWSY0NDg3w+X/Oxz+fTsWPHWqzx+/0RaxoaGpSSkhKxLhAIKBAISJJKSkqUkZFxR5vvqpibO8zNPWbnDnNrX1FfoTmO0+JrHo/nttdIUl5enkpKSlRSUqIFCxbczj7xJebmDnNzj9m5w9zcuZO5RQ2az+dTMBhsPg4Ggy1eefl8PtXX199yDQAAbSlq0DIzM1VXV6dz586pqalJ5eXlysrKiliTlZWlXbt2yXEcHT16VD169CBoAIB2FV9cXFx8qwVxcXFKT0/XunXrtG3bNo0bN05jx47V9u3bdfz4cWVmZio9PV1Hjx7V5s2bVV1drdmzZ6t3795Rn5yP9rvD3Nxhbu4xO3eYmztu5+ZxWnsDDACAewx3CgEAmEDQAAAmxHTrqzvBbbPciTa33bt3a+vWrZKkhIQEFRYWauDAgR2w084l2ty+UlNTo5dfflkvvPCCxo4d28677Hximdvhw4e1efNmhUIhJSUl6dVXX+2AnXYu0eZ25coVrV27VsFgUKFQSFOmTFFubm4H7bbz2LBhgyorK5WcnKyVK1e2eNx1F5w2FAqFnDlz5jhnzpxxvvjiC2f+/PnOf//734g1//rXv5zly5c74XDY+fe//+0sXLiwLbd0T4hlbkeOHHEuXbrkOI7jVFZWMjcntrl9ta64uNj51a9+5Xz00UcdsNPOJZa5NTY2OvPmzXM+++wzx3Ec58KFCx2x1U4llrm99957zjvvvOM4juN8/vnnzsyZM50vvviiI7bbqRw+fNg5fvy48/Of/7zVx912oU0vOXLbLHdimdvQoUOVmJgoSRoyZEjEnxXsqmKZmyR98MEHys7OVq9evTpgl51PLHPbs2ePsrOzm+8IlJyc3BFb7VRimZvH49G1a9fkOI6uXbumxMRExcXxTs+IESOaf/9qjdsutOlkW7ttVkNDQ4s1rd02qyuLZW5f9+GHH2rUqFHtsbVOLdb/3vbt26f8/Pz23l6nFcvc6urq1NjYqOLiYr300kv6xz/+0d7b7HRimdvkyZP1v//9T7Nnz9aLL76on/zkJwQtBm670KbvoTl38bZZXcntzOTQoUPasWOHli1b1tbb6vRimdvmzZs1Y8YMflP5mljmFgqFdOLECS1ZskQ3btzQ4sWLNWTIkC59r8JY5nbgwAENGDBAS5cu1dmzZ/Xaa69p2LBh6tGjR3tt857ktgttGjRum+VOLHOTpJMnT2rjxo1auHChkpKS2nOLnVIsczt+/LjWrFkjSbp48aKqqqoUFxenMWPGtOteO5NYf50mJSUpISFBCQkJGj58uE6ePNmlgxbL3Hbs2KGCggJ5PB6lp6crNTVVtbW1Gjx4cHtv957itgtt+r+p3DbLnVjmVl9frxUrVmjOnDld+jeVr4tlbuvXr2/+MXbsWBUWFnbpmEmx/zo9cuSIQqGQrl+/rpqaGvXt27eDdtw5xDI3v9+vgwcPSpIuXLig2tpapaamdsR27yluu9DmdwqprKzU7373O4XDYeXm5ur73/++tm/fLknKz8+X4zgqLS3VgQMH1K1bNxUVFSkzM7Mtt3RPiDa33/72t/r444+brzPHx8erpKSkI7fcKUSb29etX79eo0eP5mP7im1u77//vnbs2KG4uDg9+uijevzxxztyy51CtLk1NDRow4YNzR9oePLJJzV+/PiO3HKnsHr1an3yySe6dOmSkpOTNW3aNDU1NUm6sy5w6ysAgAm8Mw4AMIGgAQBMIGgAABMIGgDABIIGADCBoAEATCBoAAAT/g+FPcuV3dmPUgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig,ax=plt.subplots(1,1,figsize=(7,7))\n", "ax.set_title('landsat a17')\n", "im0=ax.imshow(np.log10(output.spectral_dens))\n", "ax.set_title('log10 of the 2-d power spectrum')\n", "divider = make_axes_locatable(ax)\n", "cax = divider.append_axes(\"bottom\", size=\"5%\", pad=0.35)\n", "out=fig.colorbar(im0,orientation='horizontal',cax=cax)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'output' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mavg_binwidth\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m5\u001b[0m \u001b[0;31m#make the kradial bins 5 pixels wide\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mannular_avg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mavg_binwidth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'output' is not defined" ] } ], "source": [ "avg_binwidth=5 #make the kradial bins 5 pixels wide\n", "output.annular_avg(avg_binwidth)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'output' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgraph_spectrum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkol_offset\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2000.\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mtitle\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Landsat {} power spectrum'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'output' is not defined" ] } ], "source": [ "output.graph_spectrum(kol_offset=2000.,title='Landsat {} power spectrum'.format(output.filename))" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# Problem -- lowpass filtering of a 2-d image\n", "\n", "For the image above, \n", "we know that the 25 meter pixels correspond to k=1/0.025 = 40 $km^{-1}$. That means that the Nyquist\n", "wavenumber is k=20 $km^{-1}$. Using that information, design a filter that removes all wavenumbers\n", "higher than 1 $km^{-1}$. \n", "\n", "1) Use that filter to zero those values in the fft, then inverse transform and\n", "plot the low-pass filtered image.\n", "\n", "2) Take the 1-d fft of the image and repeat the plot of the power spectrum to show that there is no power in wavenumbers higher than 1 $km^{-1}$.\n", "\n", "(Hint -- I used the fftshift function to put the low wavenumber cells in the center of the fft, which made it simpler to zero the outer cells. I then used ifftshift to reverse shift before inverse transforming to get the filtered\n", "image.)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "formats": "ipynb,py:percent", "notebook_metadata_filter": "all,-language_info,-toc,-latex_envs" }, "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.7.6" }, "toc": { "nav_menu": {}, "number_sections": true, "sideBar": false, "skip_h1_title": false, "toc_cell": true, "toc_position": {}, "toc_section_display": "block", "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }