gmivar_v1-sfc.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# # adapting simon's older PC-based 1dvar code to read in GMI \n",
    "# data and run oceanic retrievals. started DD 17/05/18\n",
    "#   - - working as wv retrieval 03/06/18\n",
    "#  playing with surface retrieval (jacobians first), DD late august 2018"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set up ARTS/Python environment\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "env: ARTS_INCLUDE_PATH=/home/dudavid/arts/controlfiles/\n",
      "env: ARTS_DATA_PATH=/home/dudavid/arts/arts-xml/\n",
      "env: ARTS_BUILD_PATH=/home/dudavid/arts/build/\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/dudavid/anaconda3/lib/python3.6/site-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading ARTS API from: /home/dudavid/arts/build/src/libarts_api.so\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<typhon.arts.workspace.agendas.Agenda at 0x7fbe9dc04f60>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%env ARTS_INCLUDE_PATH=/home/dudavid/arts/controlfiles/\n",
    "%env ARTS_DATA_PATH=/home/dudavid/arts/arts-xml/\n",
    "%env ARTS_BUILD_PATH=/home/dudavid/arts/build/\n",
    "\n",
    "%matplotlib inline\n",
    "import glob\n",
    "from h5py import File\n",
    "\n",
    "from typhon.arts.workspace import Workspace, arts_agenda\n",
    "ws = Workspace(verbosity=0)\n",
    "ws.execute_controlfile(\"general/general.arts\")\n",
    "ws.execute_controlfile(\"general/continua.arts\")\n",
    "ws.execute_controlfile(\"general/agendas.arts\")\n",
    "ws.execute_controlfile(\"general/planet_earth.arts\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from typhon.arts.workspace.variables import *\n",
    "\n",
    "ws.Copy( ws.abs_xsec_agenda, ws.abs_xsec_agenda__noCIA )\n",
    "ws.Copy( ws.iy_main_agenda, ws.iy_main_agenda__Emission )\n",
    "ws.Copy( ws.iy_space_agenda, ws.iy_space_agenda__CosmicBackground )\n",
    "ws.Copy( ws.propmat_clearsky_agenda, ws.propmat_clearsky_agenda__OnTheFly )\n",
    "ws.Copy( ws.ppath_agenda, ws.ppath_agenda__FollowSensorLosPath )\n",
    "ws.Copy( ws.ppath_step_agenda, ws.ppath_step_agenda__GeometricPath )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# define absorbing species and sensor (here using metmm library, used again below)\n",
    "ws.abs_speciesSet(species=[\"H2O-PWR98\",\"O2-PWR93\",\"N2-SelfContStandardType\"])\n",
    "ws.abs_lines_per_speciesSetEmpty()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ws.stokes_dim = 2 # to get V and H pol out of metmm\n",
    "ws.iy_unit = \"PlanckBT\" # equivalent: ws.StringSet( iy_unit, \"PlanckBT\" )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Atmosphere set up # (take Psfc, T from model later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[101300.  97500.  95000.  92500.  90000.  85000.  80000.  75000.  70000.\n",
      "  65000.  60000.  55000.  50000.  40000.  30000.  20000.  10000.]\n"
     ]
    }
   ],
   "source": [
    "ws.atmosphere_dim = 1  # for 1DVAR\n",
    "p = np.array([1013.,975.,950.,925.,900.,850.,800.,750.,700.,650.,600.,550.,500.,400.,300.,200.,100.])*100.0\n",
    "ws.p_grid = p[:] #0.5 * (p[1:] + p[:-1])\n",
    "print(ws.p_grid.value)\n",
    "ws.AtmRawRead( basename = \"planets/Earth/Fascod/tropical/tropical\") #tropical atmosphere assumed\n",
    "ws.AtmosphereSet1D()\n",
    "ws.AtmFieldsCalc()\n",
    "\n",
    "ws.z_surface = np.asarray(ws.z_field)[0]\n",
    "ws.t_surface = np.asarray(ws.t_field)[0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sensor definition -- GMI sensor setup w/ metmm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "ws.ArrayOfIndexCreate(\"channels\")\n",
    "# initial sensor setup, all 13 channels at low frequency EIA\n",
    "csub = [[x for x in range(9)], [x+9 for x in range(4)]] # 13 GMI channels total\n",
    "#csub = [ [3,4,5,6], [11,12] ]#  define channel indices (0-8 10-89GHz, 9-12 166-183GHz)\n",
    "nch = np.size(csub)\n",
    "#ws.channels = csub #[0,:]+csub[1,:] #list(np.ravel(csub)) # 23.8 and 183 channels only for WV\n",
    "\n",
    "ws.ArrayOfIndexCreate(\"viewing_angles\") # necessary if just using one pencil beam angle?\n",
    "#ws.viewing_angles = [1] # index of viewing angles in metmm file -- for GMI, 49.2, 52.8\n",
    "# current setup is low freqs at index 1, high freqs at index 0\n",
    "\n",
    "ws.sensor_pos  = np.array([[407e3]])#, -30.0, 180.0]]) # 407km alt, 30S, 180E\n",
    "#can take exact alt from file if desired, same with lat/lon values (does that have effect?)\n",
    "ws.sensor_time = np.array([0.0]) # same here\n",
    "ws.sensor_los  = np.array([[180.0]]) # since angle taken care of in metmm file \n",
    "# zenith angle followed by azimuth. zenith=0 is straight up, 180 is straight down. azimuth=0 is N, positive to east \n",
    "\n",
    "ws.IndexCreate(\"met_mm_accuracy\") \n",
    "ws.met_mm_accuracy = 1 # points within each bandwidth to comute abs characteristics\n",
    "\n",
    "# initial exection of metmm control files\n",
    "#ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/prepare_metmm.arts\")\n",
    "#ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_gmi.arts\") #atms.arts\")\n",
    "#ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/apply_metmm.arts\")\n",
    "#ws.execute_controlfile(\"instruments/metmm/common_metmm.arts\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform some checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ws.atmfields_checkedCalc( bad_partition_functions_ok = 1 )\n",
    "ws.abs_xsec_agenda_checkedCalc()\n",
    "ws.propmat_clearsky_agenda_checkedCalc()\n",
    "ws.atmgeom_checkedCalc()\n",
    "#ws.abs_lookupSetup() # not currently using lookup tables (on the fly instead)\n",
    "#ws.abs_lookupCalc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# define surface agendas -- unnecessary right now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#from typhon.arts.workspace import arts_agenda\n",
    "#\n",
    "#@arts_agenda\n",
    "#def surface_fastem(ws):\n",
    "#    ws.specular_losCalc()\n",
    "#    ws.InterpSurfaceFieldToPosition( out = ws.surface_skin_t, field = ws.t_surface )\n",
    "#    ws.surfaceFastem(wind_speed     = ws.wind_speed,\n",
    "#                     wind_direction = ws.wind_direction,\n",
    "#                     transmittance  = ws.transmittance)\n",
    "#    \n",
    "#    \n",
    "#ws.TessemNNReadAscii(tessem_netv, \"testdata/tessem_sav_net_V.txt\")\n",
    "#ws.TessemNNReadAscii(tessem_neth, \"testdata/tessem_sav_net_H.txt\") \n",
    "#@arts_agenda\n",
    "#def surface_tessem(ws):\n",
    "#    ws.specular_losCalc()\n",
    "#    ws.InterpSurfaceFieldToPosition( out = ws.surface_skin_t, field = ws.t_surface )\n",
    "#    ws.SurfaceTessem(wind_speed     = ws.wind_speed) # one with capital S and one without?!\n",
    "#\n",
    "## set surface agenda (i.e. either fastem or tessem) here:\n",
    "##ws.Copy(ws.surface_rtprop_agenda, surface_fastem) # choose which emis. model to run here!\n",
    "#ws.Copy(ws.surface_rtprop_agenda, surface_tessem) # choose which emis. model to run here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Designate surface vars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wind speed, SST:  9.5 287.8\n"
     ]
    }
   ],
   "source": [
    "\n",
    "ws.NumericCreate(\"wind_speed\")\n",
    "ws.NumericCreate(\"wind_direction\")\n",
    "ws.wind_speed = 9.5\n",
    "ws.surface_skin_t = 287.8\n",
    "#ws.wind_direction = 0.0 # can adjust based on analysis data later  -- default should be 0\n",
    "# -- need sensor azimuth angle, and spacecraft heading, as arts azimuth is relative to N\n",
    "#ws.salinity = .035 # default is .035 anyway\n",
    "\n",
    "# transmittance only needed for running fastem:\n",
    "#ws.VectorCreate(\"transmittance\")\n",
    "#ws.transmittance  = np.ones(ws.f_grid.value.shape) # or ones, or something else?\n",
    "print('wind speed, SST: ',ws.wind_speed.value,ws.surface_skin_t.value)#ws.t_surface.value.ravel()[0])\n",
    "\n",
    "# add surface variables for jacobian calculation -- need to keep order consistent between snames/sdata\n",
    "snames = [\"Water skin temperature\", \"Wind speed\", \"Salinity\"]\n",
    "sdata = np.array([ws.surface_skin_t.value, ws.wind_speed.value, 0.035]).reshape(3,1,1)\n",
    "#sdata = np.array([ws.t_surface.value.ravel()[0], ws.wind_speed.value, 0.035]).reshape(3,1,1)\n",
    "\n",
    "# don't copy to arts variables yet... otherwise ycalc bombs due to SurfaceDummy\n",
    "ws.Copy(ws.surface_props_names, snames)\n",
    "ws.Copy(ws.surface_props_data, sdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### want to get rid of SurfaceDummy from iy_surface agenda\n",
    "#AgendaSet( iy_surface_agenda__UseSurfaceRtprop ){\n",
    "#  SurfaceDummy\n",
    "#  iySurfaceRtpropAgenda\n",
    "#}\n",
    "\n",
    "# read tessem ascii files into arts memory\n",
    "ws.TessemNNReadAscii(tessem_netv, \"testdata/tessem_sav_net_V.txt\")\n",
    "ws.TessemNNReadAscii(tessem_neth, \"testdata/tessem_sav_net_H.txt\") \n",
    "@arts_agenda\n",
    "def iy_surface_agendaPY(ws):\n",
    "    ws.SurfaceTessem()\n",
    "    ws.iySurfaceRtpropCalc()\n",
    "    #ws.iySurfaceRtpropAgenda()\n",
    "    # don't call SurfaceDummy()!\n",
    "    \n",
    "ws.Copy(ws.iy_surface_agenda, iy_surface_agendaPY) # copy python-defined agenda to ARTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# define jacobians, run ycalc "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[169.33739751  91.60115224 199.51783459 137.41204274 236.05000015\n",
      " 220.68824268 166.40211046 271.53278258 252.60123359 283.12457811\n",
      " 283.08760084 258.52720733 272.31574126]\n",
      "[[ 4.01997415e-01 -3.31489309e-01]\n",
      " [-2.04717092e-01  9.72616606e-01]\n",
      " [ 3.56611617e-01 -2.85692054e-01]\n",
      " [-2.06388883e-01  8.65775351e-01]\n",
      " [ 3.18909736e-01 -1.82128148e-01]\n",
      " [ 3.06274073e-01 -2.96253346e-01]\n",
      " [-2.69686154e-01  9.08255492e-01]\n",
      " [ 2.80146348e-01 -1.44560026e-01]\n",
      " [-6.71921065e-02  3.80422530e-01]\n",
      " [ 3.63350921e-02 -3.15422169e-04]\n",
      " [ 1.71110613e-02  6.87726220e-04]\n",
      " [ 3.43197686e-12  8.45302559e-14]\n",
      " [ 6.54978335e-05  1.51464785e-06]]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "ws.ArrayOfIndexCreate(\"viewing_angles_1\")\n",
    "ws.viewing_angles_1 = [1] # index, defined in GMI metmm file\n",
    "ws.ArrayOfIndexCreate(\"viewing_angles_2\")\n",
    "ws.viewing_angles_2 = [0]\n",
    "\n",
    "ws.ArrayOfIndexCreate(\"channels_1\") # low freq channels EIA\n",
    "ws.channels_1 = csub[0]\n",
    "ws.ArrayOfIndexCreate(\"channels_2\") # high freq channels EIA\n",
    "ws.channels_2 = csub[1] \n",
    "\n",
    "ws.Copy(ws.viewing_angles, ws.viewing_angles_1)\n",
    "ws.Copy(ws.channels, ws.channels_1)\n",
    "\n",
    "\n",
    "# re-execution of metmm control files ( since channels subset and angle have changed)\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/prepare_metmm.arts\")\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_gmi.arts\") #atms.arts\")\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/apply_metmm.arts\") # to execute CF?\n",
    "#ws.VectorSetConstant(ws.transmittance, 1, 1.0) # needed for emis model, right size\n",
    "\n",
    "ws.jacobianInit()  # initialize jacobian quantities, then add variables\n",
    "ws.jacobianAddSurfaceQuantity(\n",
    "    g1=ws.lat_grid, g2=ws.lon_grid, quantity=snames[0])\n",
    "ws.jacobianAddSurfaceQuantity(\n",
    "    g1=ws.lat_grid, g2=ws.lon_grid, quantity=snames[1])\n",
    "#ws.jacobian_do = 1 # flag to activate clear-sky jacobian calculations, set to 1 when jacobianClose is called\n",
    "ws.jacobianClose()\n",
    "\n",
    "ws.cloudboxOff()\n",
    "ws.cloudbox_checkedCalc()\n",
    "ws.sensor_checkedCalc()\n",
    "\n",
    "ws.yCalc()  # calculate yf and jacobian matching x\n",
    "#print(ws.y.value)\n",
    "#print(ws.jacobian.value)\n",
    "\n",
    "# now run for other EIA, and append Tbs to previous EIA with yCalcAppend\n",
    "ws.Copy(ws.viewing_angles, ws.viewing_angles_2)\n",
    "ws.Copy(ws.channels, ws.channels_2)\n",
    "\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/prepare_metmm.arts\")\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_gmi.arts\") #atms.arts\")\n",
    "ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/apply_metmm.arts\") # to execute CF?\n",
    "#ws.VectorSetConstant(ws.transmittance,2,1.0) # needed for emis model, right size\n",
    "ws.yCalcAppend(jacobian_quantities_copy = ws.jacobian_quantities)\n",
    "print(ws.y.value)\n",
    "print(ws.jacobian.value)\n",
    "#\n",
    "#ws.Print(ws.jacobian_quantities,0)\n",
    "##ws.VectorAddVector( ws.yf, ws.y, ws.y_baseline )  # add baseline term\n",
    "##ws.jacobianAdjustAndTransform()\n",
    "##ws.Print(ws.jacobian,0)\n",
    "#ws.Print(ws.y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}