Geodesic vs. map projection timing for radar data

Radar projection timing test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from lmatools.coordinateSystems import RadarCoordinateSystem, MapProjection, GeographicSystem, CoordinateSystem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r = np.arange(0,10e3,5.0, dtype='float')[:,None,None]\n",
    "az = np.arange(0, 360, 1, dtype=float)[None,:,None]\n",
    "el = np.arange(0,30,2,dtype=float)[None,None,:]\n",
    "\n",
    "r, az, el = np.meshgrid(r, az, el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(360, 2000, 15) (360, 2000, 15) (360, 2000, 15)\n",
      "10800000\n"
     ]
    }
   ],
   "source": [
    "print r.shape, az.shape, el.shape\n",
    "print r.shape[0] * az.shape[1] * el.shape[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ctrlat, ctrlon, ctralt = 33.0, -101.0, 0.0\n",
    "radar = RadarCoordinateSystem(ctrlat, ctrlon, ctralt)\n",
    "mproj = MapProjection(projection='aeqd', ctrLat=ctrlat, ctrLon=ctrlon, lat_ts=ctrlat, \n",
    "                            lon_0=ctrlon, lat_0=ctrlat, lat_1=ctrlat)\n",
    "geo = GeographicSystem()\n",
    "\n",
    "lon, lat, alt = radar.toLonLatAlt(r,az,el)\n",
    "X, Y, Z = geo.toECEF(lon,lat,alt)\n",
    "mx, my, mz = mproj.fromECEF(X, Y, Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 5.51 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "lon, lat, alt = radar.toLonLatAlt(r,az,el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.71 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "X, Y, Z = geo.toECEF(lon,lat,alt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 5.64 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "mx, my, mz = mproj.fromECEF(X, Y, Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pyproj as proj4\n",
    "class MapProjectionSimple(CoordinateSystem):\n",
    "    \"\"\"Map projection coordinate system. Wraps proj4, and uses its projecion names. Defaults to \n",
    "        equidistant cylindrical projection\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, projection='eqc', ctrLat=None, ctrLon=None, ellipse='WGS84', datum='WGS84', **kwargs):\n",
    "        self.projection = proj4.Proj(proj=projection, ellps=ellipse, datum=datum, **kwargs)\n",
    "        self.ctrLat=ctrLat\n",
    "        self.ctrLon=ctrLon\n",
    "        self.ctrAlt=0.0\n",
    "        self.geoCS = GeographicSystem()\n",
    "\n",
    "    def fromLonLatAlt(self, lon,lat,alt):\n",
    "        projectedData = np.array(proj4.transform(self.geoCS.WGS84lla, self.projection, lon, lat, alt))\n",
    "        if len(projectedData.shape) == 1:\n",
    "            px, py, pz = projectedData[0], projectedData[1], projectedData[2]\n",
    "        else:\n",
    "            px, py, pz = projectedData[0,:], projectedData[1,:], projectedData[2,:]\n",
    "        return px, py, pz\n",
    "\n",
    "\n",
    "projsimple = MapProjectionSimple(projection='aeqd', ctrLat=ctrlat, ctrLon=ctrlon, lat_ts=ctrlat, \n",
    "                            lon_0=ctrlon, lat_0=ctrlat, lat_1=ctrlat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 3.39 s per loop\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "mx, my, mz = projsimple.fromLonLatAlt(lon, lat, alt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}