am2var_v1-sfcr.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\n",
    "\n",
    "# revisited DD mid Dec 2018 for running 1dvar sfc retrieval\n",
    "# trying retrieval with AMSR2 metmm sensor setup 17/12/18"
   ]
  },
  {
   "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 0x7f1d85c0c588>"
      ]
     },
     "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
   },
   "outputs": [],
   "source": [
    "#import typhon\n",
    "#typhon.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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 -- AMSR2 sensor setup w/ metmm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# channels using:  8\n"
     ]
    }
   ],
   "source": [
    "ws.ArrayOfIndexCreate(\"channels\")\n",
    "# initial sensor setup, all 13 channels at low frequency EIA\n",
    "nlo = 14\n",
    "#csub = [[x for x in range(nlo)], [x+nlo for x in range(2)]] # 16 AMSR2 channels total\n",
    "csub = [[0,1,2,3,4,5,6,7],[]] # so 6.9V/H, 7.3V/H, 10V/H, 18V/H\n",
    "#csub = [ [3,4,5,6], [11,12] ]#  define channel indices (0-8 10-89GHz, 9-12 166-183GHz)\n",
    "nch = np.size(np.ravel(csub[0]+csub[1]))\n",
    "print('# channels using: ',nch)\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 AMSR2 only B-scan is different\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 = 2 # points within each bandwidth to comute abs characteristics\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perform some checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "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": 10,
   "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": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "wind speed, SST:  3.5 290.2\n"
     ]
    }
   ],
   "source": [
    "\n",
    "ws.NumericCreate(\"wind_speed\")\n",
    "ws.NumericCreate(\"wind_direction\")\n",
    "ws.wind_speed = 3.5\n",
    "ws.surface_skin_t = 290.2\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",
    "#print(ws.f_grid.value)\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\", \"Wind direction\"]\n",
    "sdata = np.array([ws.surface_skin_t.value, ws.wind_speed.value, 0.035, ws.wind_direction.value]).reshape(4,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 [old issue]\n",
    "ws.Copy(ws.surface_props_names, snames)\n",
    "ws.Copy(ws.surface_props_data, sdata)\n",
    "\n",
    "\n",
    "nrvar = 2 # num retrieved variables\n",
    "rvars = np.zeros(nrvar)\n",
    "#rvars[:] = [0.8, 2.2] # corresponding to variable order set above\n",
    "rvars[:] = [4.2, 3.9] # corresponding to variable order set above\n",
    "\n",
    "# set Xa (formally set below)\n",
    "xa = np.zeros(nrvar) # if supplying user-defined a priori vector\n",
    "xa[:] = [ws.surface_skin_t.value, ws.wind_speed.value] # just perturbations from values defined above\n",
    "#xa[:] = [ws.surface_skin_t.value + .42, ws.wind_speed.value - 0.78] # just perturbations from values defined above"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# define jacobians, run ycalc "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<typhon.arts.workspace.agendas.Agenda at 0x7f1d85c1d320>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "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",
    "ws.Copy(ws.viewing_angles, ws.viewing_angles_1)\n",
    "ws.Copy(ws.channels, ws.channels_1)\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_amsr2.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",
    "# NO NEED TO RUN JACOBIAN INIT IF USING RETRIEVAL DEF INIT\n",
    "#ws.jacobianInit()  # initialize jacobian quantities, then add variables\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[6.8375e+09 7.0125e+09 7.2125e+09 7.3875e+09 1.0625e+10 1.0675e+10\n",
      " 1.8650e+10 1.8750e+10]\n"
     ]
    }
   ],
   "source": [
    "#### want to get rid of SurfaceDummy from iy_surface agenda\n",
    "#AgendaSet( iy_surface_agenda__UseSurfaceRtprop ){\n",
    "#  SurfaceDummy\n",
    "#  iySurfaceRtpropAgenda\n",
    "#}\n",
    "\n",
    "# if using tessem... 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",
    "\n",
    "# if using surfaceFastem, need to provide transmittances and wind direction:\n",
    "ws.VectorCreate( \"transmittance\" )\n",
    "##nelemGet( v = f_grid )\n",
    "print(ws.f_grid.value)\n",
    "ws.transmittance = np.ones( ws.f_grid.value.shape )\n",
    "#ws.Copy(ws.tranny, tran)\n",
    "#ws.VectorSetConstant( ws.transmission_field, ws.f_grid.value.shape, 0.9 )\n",
    "##and call the method as\n",
    "#SurfaceFastem( transmittance = tr )\n",
    "\n",
    "\n",
    "@arts_agenda\n",
    "def iy_surface_agendaPY(ws):\n",
    "    ws.SurfaceFastem( transmittance = ws.transmittance )   #### DOES THIS WORK\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": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[157.65357936  83.79344135 158.13948819  84.30832419 163.10230934\n",
      "  89.89433944 193.41729417 132.96342163]\n",
      "[[ 0.74461149 -0.24161547]\n",
      " [ 0.15633142  0.54459729]\n",
      " [ 0.73693003 -0.24338306]\n",
      " [ 0.15401885  0.55124132]\n",
      " [ 0.63389294 -0.25928305]\n",
      " [ 0.12019723  0.59868726]\n",
      " [ 0.34929299 -0.25640154]\n",
      " [ 0.02391223  0.64347634]]\n"
     ]
    }
   ],
   "source": [
    "#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.jacobianClose()\n",
    "ws.retrievalDefInit()  # initialize jacobian quantities, then add variables\n",
    "\n",
    "sx_cov = np.diag(rvars) # CHANGE TO MATCH SIZE OF Xa!!\n",
    "\n",
    "ws.retrievalAddSurfaceQuantity(\n",
    "    g1=ws.lat_grid, g2=ws.lon_grid, quantity=snames[0])\n",
    "ws.covmat_sxAddBlock(block = sx_cov[:1, :1])\n",
    "\n",
    "ws.retrievalAddSurfaceQuantity(\n",
    "    g1=ws.lat_grid, g2=ws.lon_grid, quantity=snames[1])\n",
    "ws.covmat_sxAddBlock(block = sx_cov[1:, 1:])\n",
    "\n",
    "ws.Copy(ws.xa, xa)\n",
    "# print Xa to terminal:\n",
    "ws.Print(ws.xa,0) # Print Xa\n",
    "\n",
    "## set channel errors, currently just constant across channels\n",
    "chan_err = 0.6 # for now all the same, in K\n",
    "\n",
    "se_cov = np.diag(np.zeros(nch) + chan_err**2 ) # set channel assumed errors\n",
    "#print('Se diag ',np.diag(se_cov))\n",
    "ws.covmat_seSet(se_cov)\n",
    "\n",
    "ws.retrievalDefClose()\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",
    "am_tb = ws.y.value\n",
    "print(ws.jacobian.value)\n",
    "\n",
    "# now run for other EIA, and append Tbs to previous EIA with yCalcAppend\n",
    "if csub[1] != []: \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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# setup for retrieval and sample jacobian seem right, time to run a sample retrieval:\n",
    "\n",
    "@arts_agenda\n",
    "def inversion_iterate_agenda(ws):\n",
    "    ws.Ignore(ws.inversion_iteration_counter) # on simon's advice\n",
    "    ws.x2artsAtmAndSurf()  # map x to ARTS's variables\n",
    "    \n",
    "    # NOTE: IF WANTING TO CHANGE CHANNELS/ANGLES, NEED TO USE arts COMMANDS ONLY!\n",
    "    \n",
    "    ws.Copy(ws.viewing_angles, ws.viewing_angles_1)\n",
    "    ws.Copy(ws.channels, ws.channels_1)\n",
    "    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/prepare_metmm.arts\")\n",
    "    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_amsr2.arts\") \n",
    "    #ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_gmi.arts\") \n",
    "    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/apply_metmm.arts\") \n",
    "    ws.yCalc() \n",
    "\n",
    "    #if csub[1] != []: \n",
    "    #    # aka if any 166/183 channels exist in channel selection, run other angle\n",
    "    #    ws.Copy(ws.viewing_angles, ws.viewing_angles_2)\n",
    "    #    ws.Copy(ws.channels, ws.channels_2)\n",
    "    #    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/prepare_metmm.arts\")\n",
    "    #    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/sensor_gmi.arts\")\n",
    "    #    ws.execute_controlfile(\"instruments/metmm/sensor_descriptions/apply_metmm.arts\")\n",
    "    #    ws.yCalcAppend(jacobian_quantities_copy = ws.jacobian_quantities)\n",
    "        \n",
    "    ws.Copy(ws.y_baseline,np.zeros(nch)) # kinda silly but necessary?\n",
    "    ws.VectorAddVector( ws.yf, ws.y, ws.y_baseline )  # add baseline term\n",
    "    \n",
    "    ws.jacobianAdjustAndTransform()\n",
    "    # last line takes care of some fixes needed to get the jacobian right for iterative solutions\n",
    "\n",
    "ws.Copy(ws.inversion_iterate_agenda, inversion_iterate_agenda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ws.inversion_iterate_agenda.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in GMI TBs to use as observation vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Import GMI L1C file with brightness temperatures (TBs)\n",
    "#paf = '/home/dudavid/Dendrite/Dendrite/SatData/GMI/L1C/'\n",
    "##l1c = glob.glob(paf+'*.HDF5')\n",
    "##print(l1c[0])\n",
    "##\n",
    "##for f in l1c[0:1]:\n",
    "##    gfi = File(f)\n",
    "##    tb1 = gfi['S1/Tc'][:]   # given as [scans, pix], ie 2962x221\n",
    "##    tb2 = gfi['S2/Tc'][:]\n",
    "##    #print(np.shape(tbs),type(tbs))\n",
    "##    tbs = np.append(tb1, tb2, axis=2) #gfi['S2/Tc'][:],axis=1)\n",
    "##    print(np.shape(tbs))\n",
    "#    #lat = gfi['S1/Latitude'][:]\n",
    "#    #lon = gfi['S1/Longitude'][:]\n",
    "#    #eia = np.append(gfi['S1/incidenceAngle'][:],gfi['S2/incidenceAngle'][:],axis=2)\n",
    "#    #print(np.shape(eia)) # EIA is one for S1 channels, one for S2 (approx 52.8, 49.2 deg)\n",
    "#    # READ THESE IN BUT NOT CURRENTLY USED\n",
    "#    # can also import SC data, time data, etc from L1C file\n",
    "#\n",
    "#gfi = File(paf+'1C-R.GPM.GMI.XCAL2016-C.20180418-S160617-E173849.023507.V05A.HDF5') # hand-picked orbit\n",
    "#tb1 = gfi['S1/Tc'][:]\n",
    "#tb2 = gfi['S2/Tc'][:]\n",
    "#tbs = np.append(tb1, tb2, axis=2) \n",
    "#print(np.shape(tbs))\n",
    "#lat = gfi['S1/Latitude'][:]\n",
    "#lon = gfi['S1/Longitude'][:]\n",
    "#sel1 = np.logical_and(lon>-130, lon<-110)\n",
    "#sc,px = np.indices([lat.shape[0],lat.shape[1]])\n",
    "#sel2 = np.logical_and(lat[sel1]>0, lat[sel1]<20)\n",
    "## scans where we're in about the right lat/lon range specified above:\n",
    "#scrange = np.round([np.mean(sc[sel1][sel2])-20, np.mean(sc[sel1][sel2])+20])\n",
    "#gmi_tbs = tbs[int(scrange[0]),90,:][:] # sure\n",
    "#gmi_tb = gmi_tbs[np.ravel(csub[0]+csub[1])]\n",
    "#\n",
    "## ...or use something modified by hand to be like the forward modeled one given the prior...\n",
    "##gmi_tb = np.array([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",
    "#print(gmi_tb)\n",
    "##gmi_tb -= .71\n",
    "##gmi_tb[0] += .4\n",
    "##gmi_tb[2] += .7\n",
    "\n",
    "ws.Copy(ws.y, am_tb)\n",
    "#am_tb[0]+=0.8\n",
    "am_tb[1]+=1.1\n",
    "am_tb[7]+=2.9\n",
    "#ws.Copy(ws.y, gmi_tb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ARTS[94475840951952]:     5        1.96566       0.188402        1.77725     0.00901448\n",
      "ARTS[94475840951952]: --------------------------------------------------------------------------------\n",
      "ARTS[94475840951952]: \n",
      "ARTS[94475840951952]: Total number of steps:            5\n",
      "ARTS[94475840951952]: Final scaled cost function value: 1.96566\n",
      "ARTS[94475840951952]: OEM computation converged.\n",
      "ARTS[94475840951952]: \n",
      "ARTS[94475840951952]: Elapsed Time for Retrieval:                       0.0265587\n",
      "ARTS[94475840951952]: Time in inversion_iterate Agenda (No Jacobian):   0.00367871\n",
      "ARTS[94475840951952]: Time in inversion_iterate Agenda (With Jacobian): 0.0226066\n",
      "ARTS[94475840951952]: \n",
      "ARTS[94475840951952]:                                       ----                                      \n",
      "ARTS[94475840951952]: \n",
      "\n"
     ]
    }
   ],
   "source": [
    "ws.OEM(method=\"gn\", #lm\",\"gn\"\n",
    "    max_iter=100,\n",
    "    display_progress=1,\n",
    "    max_start_cost=1e5,\n",
    "    lm_ga_settings=np.array([10.0,2.0,3.0,10000000000.0,1.0,1.0]))\n",
    "#ws.x2artsStandard() # convert from ARTS coords back to user-defined grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ws.oem_errors.print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1d84145d30>]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl41OW5//H3nT1kJckkgbBDFkgC\nCBHcAXEBREDAU1t72np66lm6/uq+K4JiXdraerRWPbWnra1gUdQEUQRExSWoIWHJQiAQyEZCQkIS\nsj2/PzJQigkJZGa+s9yv65orM5Nv5rnNhZ/55p7n+zxijEEppZRv8bO6AKWUUq6n4a+UUj5Iw18p\npXyQhr9SSvkgDX+llPJBGv5KKeWDNPyVUsoHafgrpZQP0vBXSikfFGB1Ab2Ji4szo0aNsroMpZTy\nKNu2bTtsjLH1dZzbhv+oUaPIzc21ugyllPIoIlLWn+O07aOUUj5Iw18ppXyQhr9SSvkgDX+llPJB\nGv5KKeWDNPyVUsoHafgrpZQP0vBXSik38tq2cl7NPeD0cTT8lVLKTRhjeGZjCW98ddDpY2n4K6WU\nmyiqaqL08DHmZgxx+lga/kop5Say8ysQgavTE50+loa/Ukq5iZyCCqaNisEWEez0sTT8lVLKDZRU\nN1JU1cTcDOef9YOGv1JKuYWc/EoA5rig3w8a/kop5RayCyqZOnIwiVEhLhlPw18ppSy27/AxdlUc\ndVnLBzT8lVLKcjkF3S2fuZmuafmAhr9SSlkup6CCScOjSYoOddmYGv5KKWWhA3XNbC9vcGnLBzT8\nlVLKUutOtHw0/JVSyndkF1SQPjSSkbFhLh1Xw18ppSxS0dDCl/vrmefCD3pP0PBXSimLWNXyAQ1/\npZSyTE5+JWmJEYyxhbt8bA1/pZSyQPXRVj4vq3PJ8s090fBXSikLvLOjEmNgbqbrWz6g4a+UUpbI\nzq9krC2M5HjXt3xAw18ppVyutuk4n+6tZV7mEETEkho0/JVSysXW76yiy2BZvx8cFP4i8pKIVItI\nQS/fnykiDSLylf12vyPGVUopT5SdX8Go2EGMHxJhWQ2OOvP/AzCnj2O2GGMm22/LHDSuUkp5lCPH\n2vh4Ty1zLWz5gIPC3xjzAVDniNdSSilv9u6uKjq7DPMsbPmAa3v+F4pInojkiEh6TweIyM0ikisi\nuTU1NS4sTSmlXCMnv4Jhg0PJSIq0tA5Xhf8XwEhjzCTgN8DrPR1kjHneGJNljMmy2WwuKk0ppVyj\noaWdD0sOMzcj0dKWD7go/I0xR40xTfb72UCgiMS5YmyllHIXG3ZV0d5pXLpjV29cEv4ikij2tzkR\nmWYft9YVYyullLvIKahkSFQIk4dFW10KAY54ERF5BZgJxIlIOfAAEAhgjHkOWAr8l4h0AC3ADcYY\n44ixlVLKEzQd72BzUQ03Th+Bn5+1LR9wUPgbY77Zx/d/C/zWEWMppZQnen93NW0dXZas3d8TvcJX\nKaVcICe/AltEMFNHDLa6FEDDXymlnK65rYONhdXMSU90i5YPaPgrpZTTbSqsobW9y7Llm3ui4a+U\nUk6WnV9BbFgQ00bFWF3KSRr+SinlRK3tnWzcXc1V6YkE+LtP5LpPJUop5YU+KKrhWFsn89yo5QMa\n/kop5VQ5BZVEDwrkgjGxVpfyTzT8lVLKSY53dPLeziquHJ9AoBu1fEDDXymlnOajksM0Hu9wmwu7\nTqXhr5RSTpKdX0lESAAXjXOvlg9o+CullFO0d3bxrr3lExzgb3U5X6Phr5RSTrB1Ty0NLe1usXxz\nTzT8lVLKCXIKKggL8ufSZPfcukTDXymlHKyjs4t3dlQxe3wCIYHu1/IBDX+llHK4z/bWUXesjbkZ\n7nVh16k0/JVSysGyCyoIDfRnZmq81aX0SsNfKaUcqLPLsK6gillpNkKD3LPlAxr+SinlUNvKjnC4\n6ThzM9xzls8JGv5KKeVA2fkVBAf4MSvNfVs+oOGvlFIO09VlWFdQyYwUG+HBDtki3Wk0/JVSykG+\nPFBP5dFWt9qxqzca/kop5SA5+RUE+guzxydYXUqfNPyVUsoBjDHkFFRyabKNyJBAq8vpk4a/Uko5\nwPbyBg7Wt7j1hV2n0vBXSikHyCmoJMBPuHKC+7d8QMNfKaUGrLvlU8FF4+KIHhRkdTn94pDwF5GX\nRKRaRAp6+b6IyNMiUiIi20VkiiPGVUopd7Cz4ihltc3M85CWDzjuzP8PwJwzfH8ukGy/3Qw866Bx\nlVLKcjn5lfgJHtPyAQeFvzHmA6DuDIcsBP5oun0CRIuIe1/7rJRS/WCMITu/ggvGxBIbHmx1Of3m\nqp5/EnDglMfl9ueUUsqjFVU1UXr4mNvu2NUbV4W/9PCc+dpBIjeLSK6I5NbU1LigLKWUGpjs/ApE\n4Op0z2n5gOvCvxwYfsrjYcCh0w8yxjxvjMkyxmTZbDYXlaaUUuduXUEl54+KIT4ixOpSzoqrwn8t\n8B37rJ8LgAZjTIWLxlZKKacoqW6isKrRo2b5nOCQZedE5BVgJhAnIuXAA0AggDHmOSAbmAeUAM3A\nTY4YVymlrLSuoPscdo6br93fE4eEvzHmm3183wA/dMRYSinlLrLzK5kyIprEKM9q+YBe4auUUudk\n3+Fj7Kw4yjwPm+Vzgoa/Ukqdg5yCSgDmeGC/HzT8lVLqnKwrqGDSsCiGDR5kdSnnRMNfKaXOUvmR\nZvLKGzzuwq5TafgrpdRZWmdv+XjK2v090fBXSqmzlJ1fwYQhkYyMDbO6lHOm4a+UUmehoqGFL/bX\nM88DNmk/Ew1/pZQ6CydbPh7c7wcNf6WUOis5+ZWkJkQw1hZudSkDouGvlFL9VN3Yyudldcz18JYP\naPgrpVS/vbOjCmPw2Kt6T6Xhr5RS/ZSTX8FYWxjJ8Z7d8gENf6WU6pfapuN8UlrL3IwhiPS0P5Vn\n0fBXSql+WL+zii6DV/T7QcNfKaX6JTu/gpGxg5gwJNLqUhxCw18ppfpQ39zG1j3e0/IBDX+llOrT\nuzur6OgyHn9V76k0/JVSqg85BZUkRYeSmRRldSkOo+GvlFJncLS1nS3FNczLTPSalg9o+Cul1Blt\n2FVFe6fxyE3az0TDXymlziA7v5LEyBDOGx5tdSkOpeGvlFK9aDreweaiGuZkJOLn5z0tH9DwV0qp\nXr2/u5q2ji6vWMvndBr+SinVi3UFFdgigpk6crDVpTichr9SSvWgua2DjbtrmJOeiL+XtXxAw18p\npXq0ubCGlvZOj96k/Uw0/JVSqgfZBZXEhAUxbXSM1aU4hUPCX0TmiEihiJSIyJ09fP97IlIjIl/Z\nb//uiHGVUsoZWts7eX9XFVenJxDg753nyAEDfQER8QeeAa4EyoHPRWStMWbnaYf+zRjzo4GOp5RS\nzvZBUQ3H2jqZ62UXdp3KEW9p04ASY0ypMaYN+Cuw0AGvq5RSllhXUElUaCAXjo21uhSncUT4JwEH\nTnlcbn/udEtEZLuIrBaR4T29kIjcLCK5IpJbU1PjgNKUUursHO/o5N1dVVw1IYFAL235gGPCv6c5\nUOa0x28Co4wxE4H3gJd7eiFjzPPGmCxjTJbNZnNAaUopdXY+LqmlsbXDKy/sOpUjwr8cOPVMfhhw\n6NQDjDG1xpjj9oe/B6Y6YFyllHK47PwKIoIDuGic97Z8wDHh/zmQLCKjRSQIuAFYe+oBInLqW+gC\nYJcDxlVKKYdq7+xi/c4qrpiQQHCAv9XlONWAZ/sYYzpE5EfAO4A/8JIxZoeILANyjTFrgZ+IyAKg\nA6gDvjfQcZVSytG27qmloaXday/sOtWAwx/AGJMNZJ/23P2n3L8LuMsRYymllLPkFFQQFuTPZSne\n/5mj936UrZRSZ6Gjs4v1O6q4fHwCIYHe3fIBDX+llALgs3111B5rY54PtHxAw18ppQDIya8kNNCf\nmanxVpfiEhr+Simf19llWLejkpmpNkKDvL/lAxr+SinFtrIj1DQeZ66XX9h1Kg1/pZTPy86vICjA\nj8vTfKPlAxr+Sikf19VlWFdQyYwUG+HBDpn97hE0/JVSPu2r8noqj7YyL9M3ZvmcoOGvlPJpOfkV\nBPoLs8cnWF2KS2n4K6V8ljGG7PxKLhkXR2RIoNXluJSGv1LKZ+UfbOBgfYtPzfI5wXc+3XAzR1vb\n+bjkMJuLamhoaefSZBszU20MiQq1ujSlfEZ2fiUBfsJVE3yr5QMa/i7T1WXYWXGUzUU1bC6sYdv+\nI3R2GSKCAwgPCSA7vxKAtMQIZqbGMyvVxpSRg716JyGlrGSMIaegggvHxhI9KMjqclxOw9+J6o61\nsaW4O+w/KK7hcFMbABlJkfznjDHMSInnvBHRBPgJxdVNbNxdzabCGl7YUspzm/cQERLApclxzEyN\nZ2aKjfjIEIv/i5TyHrsqGimrbeY/Z4y1uhRLaPg7UEdnF3nlDWwurGZzUQ3bDzZgDAweFMhlKTZm\npNi4NNmGLSL4az+bkhBBSkIE/zFjLI2t7XxUcphNhTVsLKw++VdBRlIks1LjmZlqY/Lwwfj79bSD\nplKqP3IKKvATfLLlAyDGnL7drnvIysoyubm5VpfRp6qjrd2tnKIaPiw+TENLO34Ck4dHMyMlnhmp\nNjKTos45qI0x7KpoZGNh9T+1i6IHBXKZ/XOCGSk2YsO//oailOqZMYbZT20mMTKEv/zgAqvLcSgR\n2WaMyerrOD3zP0ttHV3kltWd7N3vrmwEID4imKsmJDAj1cYl4+Ic1kMUESYMjWTC0Eh+OGscDc3t\nbCmpYePu7jectXmHEIGJw6KZlWpjVmo8mUlR+OlfBUr1qri6idKaY9x00SirS7GMhn8/HKhrZnNR\nDZsKa9i65zDH2joJ9BeyRsZw59w0ZqTYSEuMQMT5gRs1KJD5E4cyf+JQuroMOw4dZWNhNRsLq/n1\nhmJ+9V4xsWFBzEixMTMtnsuSHfdGpJS3yM6vQASuTvetq3pPpeHfg9b2Tj4prT3ZzimtOQZAUnQo\ni85LYkaKjYvGxVm+Doifn5A5LIrMYVH8ZHbyyQ+YN+7ufjP4+5cH8RM4b8RgZqXamJkaT/rQSJe8\nSSnlznLyKzl/ZIxPT6LQ8Ke7/7en5tjJsP+0tJbjHV0EB/hxwZhYvj19JDNSbYyJC3Pr4IwJC2Lh\n5CQWTk6is8uQV17PpsIaNhVW88T6Ip5YX0R8RDAzUmzMSovnkmTfu6pRqT01TRRWNfLAtROsLsVS\nPhv+Tcc7+LjkMJvsvfuD9S0AjLGFcaM97KePjvHYvTz9/YQpIwYzZcRgfn5lCjWNx/mgqHv20Ds7\nKlm1rZwAP2HqyMHd1xWk2UhNcE3rSikrrSvonj03x0e2a+yNz8z2OTFrpvvsvprcfUfo6DKEBflz\n0bg4ZtinYg6PGeSwMd1VR2cXXx6oZ1NhNRt317Cz4igAQ6JCuq8pSLVxsRu0tZRyhmue3kJwgB9/\n/++LrS7FKXS2D1Df3MaW4u4lFD4oqqG68TgA44dE8u+XjmFGio2pIwcTFOBbV9EG+Ptx/qgYzh8V\nw21Xp1HZ0Mrmou4LzN7MO8Qrn+0n0F+YNjrm5HUFY23h+leB8nhltcfYcego914z3upSLOd14X/k\nWBt/3FrGpqJq8g7U02UgKjSQS5LjmJli47IUGwk+/CFPTxKjQvjG+SP4xvkjaOvoYlvZETYVdr8Z\nLH97F8vf3sWwwaHMSo3n8rR4ZqTYdCqp8kg59paPL8/yOcHr2j4NLe1kLX+XCUOjTrZyJg2LIkDX\nyDknB+tbTraHPt5zmOa2Tn5y+Th+flWq1aUpddYW/vZDDLD2R5dYXYrT+GzbJyo0kC/uu5IIncXi\nEEnRodw4fSQ3Th/J8Y5O7li9nWc27eGq9EQykqKsLs+rlNY0MdrNZ5R5svIjzeSVN3DHnDSrS3EL\nDjkdFpE5IlIoIiUicmcP3w8Wkb/Zv/+piIxyxLi90eB3juAAfx5ckM7gQUHcvno77Z1dVpfkNdbm\nHeLyJzfzm/dLrC7Fa52Y5TPXx2f5nDDg8BcRf+AZYC4wAfimiJw+gfb7wBFjzDjgl8BjAx1XWSN6\nUBDLF2Wws+Iov9u8x+pyvEJ1Yyv3v1GAv5/w240l7Dt8zOqSvFJOQSUThkQyKi7M6lLcgiPO/KcB\nJcaYUmNMG/BXYOFpxywEXrbfXw3MFv3b1mPNyUhk/sQhPL2hhKKqRqvL8WjGGO56LZ+Wtk7+9P3p\nBPv7cd8bBbjrZ3GeqrKhlW1lR3xuk/YzcUT4JwEHTnlcbn+ux2OMMR1AAxDrgLGVRR5akE54SAC3\nrd5OZ5cG1blava2cDburue3qVC4cG8ttc1LZUnyYt7ZXWF2aV1lX0P37nJPhe9s19sYR4d/TGfzp\nadCfYxCRm0UkV0Rya2pqHFCacpbY8GAeXJBO3oF6Xvpwr9XleKRD9S0se3Mn00bF8G8Xjwbgxukj\nmTgsimVv7eRoa7vFFXqP7IJKUhLCGRcfbnUpbsMR4V8ODD/l8TDgUG/HiEgAEAXUnf5CxpjnjTFZ\nxpgsm83mgNKUM107cQhXTkjgifWF7NU+9VkxxnDHa9vp6DI8fv3Ek9dN+PsJKxZlUtt0nCffKbS4\nSu9wqL6Fz/fVMVfP+v+JI8L/cyBZREaLSBBwA7D2tGPWAt+1318KvG+0qenxRITlizIIDvDjjtXb\n6dL2T7/95bP9bCk+zN3z0hgZ+88fQGYOi+I7F47ij5+Usb283qIKvYMxhvteLyDI348lU4ZZXY5b\nGXD423v4PwLeAXYBrxpjdojIMhFZYD/sRSBWREqAnwNfmw6qPFNCZAj3zZ/AZ/vq+NOnZVaX4xEO\n1DWz4u1dXDwulhunj+zxmFuuSsEWHszda/L1M5UBeO2Lgyc/UxkR6/3rdp0Nh8zzN8ZkG2NSjDFj\njTEr7M/db4xZa7/faoy53hgzzhgzzRhT6ohxlXtYOnUYl6XYWJmzmwN1zVaX49a6ugy3rsrDT4Rf\nLJ3U6zIZESGBPHBtOgUHj/J/W/e5tEZvUdnQykNv7uD8UYO5yf6ZivoHXfNADZiI8OjiTAS46+/5\nOk3xDF7euo9P99Zx//wJJEWHnvHYeZmJzEix8cT6IqqOtrqmQC9hjOGuv3dfiPj40knnvIe2N9Pw\nVw6RFB3KXfPG82HJYV7NPdD3D/ig0pomHlu3m1mpNq7P6rv/LCIsW5hOe2cXy97a6YIKvceqbeVs\nLKzhjjlpelFXLzT8lcN8a9oILhgTw/K3dlHZoGeqp+q0t3uCA/xZuWRiv9fvGRkbxo8vH8fb2yvY\nVFjt5Cq9Q0VDCw+/uZNpo2P47oWjrC7HbWn4K4fx8xMeWzKR9q4u7lmj7Z9T/X5LKV/sr+ehBeln\nvaT4Dy4bw1hbGPe/sYPW9k4nVegduqfQ5ndPoV06UZcePwMNf+VQI2PDuO3qNDbsruaNr06/3MM3\nFVU18tT6Iq5OT2Dh5KFn/fPBAf4sX5TJ/rpmntmoC7+dyau5B/igqIY75359Cq36Zxr+yuG+d9Eo\npoyI5sE3d1Dd6Nvtn/bOLm55NY/wkABWXJd5zss1Xzg2lsVTknhu8x5KqpscXKV3OFjfwvK3dnHB\nmBj+9YKep9Cqf9DwVw7n79c9jbG5rZMH3thhdTmWenbTHvIPNrB8UQZx4cEDeq27541nUFAA976u\nLbXTGWO487XtdBrDL5b0PoVW/YOGv3KKcfHh/L8rUsgpqCQ73zcXKdtxqIGnNxSzYNJQ5mUOfGmB\nuPBg7pybxieldaz58qADKvQef/38AFuKD3PX3DS9mKufNPyV0/zg0tFkJkVx/xsF1B1rs7oclzre\n0cktr+YxOCyIZQvTHfa638gazpQR0ax4exf1zb71O+1N+ZFmlr+1k4vG9n7FtPo6DX/lNAH+fjx+\n/UQaWtpZ9qZvtX+e3lDM7spGVi7OJHpQkMNe189PWHFdJvUt7Ty2Thd+O7FAHsBjS3R2z9nQ8FdO\nlZYYyQ9njeP1rw7x3s4qq8txia8O1PPspj0snTqM2eMTHP7644dE8v1LRvPKZ/vZVva1xXF9yp8/\n3c9HJbXcfc14hsdou+dsaPgrp/vvmeNIS4zgntfzaWjx7jXqW9s7ueXVr0iIDOH+a0/fzdRxfjo7\nmaFRIdyzpsBn91I+UNfMo9m7uGRcHN+aNsLqcjyOhr9yuqAAPx5fOonDTW088vYuq8txqifXF7Kn\n5hiPLZlIZEig08YJCw7gwQXp7K5s5A8f7XPaOO6qq6u73SMirFxy7lNofZmGv3KJzGFR3HzZGP6W\ne4Atxd65S9vn++p44cO93Dh9BJelOH8zoqvSE7lifAK/fK+Ig/UtTh/Pnfz50zI+3lPLPdeMZ9hg\nbfecCw1/5TI/nZ3MGFsYd76WT9PxDqvLcajmtg5uXZXHsMGh3D1vvMvGfXDBBIyBh9b6zgfq+2ub\neTRnN5cmx3HD+cP7/gHVIw1/5TIhgf48vnQihxpa+MW63VaX41Arc3ZTVtvM40snERYc4LJxhw0e\nxM+uSGb9zire9YEP1Lu6DLetzsNfuteR0nbPudPwVy41dWQMN100mj9uLePT0lqry3GIj0oO88et\nZdx08SguGBPr8vH/7ZLRpCZE8ODaHTS3eddfVKf7v0/K+HRvHffOH8/QPvZDUGem4a9c7tarUxgR\nM4g7XttOS5tnr1LZ2NrO7au3MzoujNuvTrOkhkB/P1Zcl8HB+hZ+vaHYkhpcoaz2GCtzdjMjxca/\nZGm7Z6A0/JXLDQoK4LElE9lX28xT73r2hUor3t5FRUMLT1w/idAgf8vqyBoVww3nD+fFLXvZXXnU\nsjqcpbvds50Af53d4yga/soSF46N5dsXjODFD/fyxf4jVpdzTjYWVvPXzw/wg8vGMHXkYKvL4Y45\naUSGBnLvmgK6vGzT95e37uOzvXXcN38CQ6K03eMIGv7KMnfOHc+QqFBuX72d4x2e1f5paG7nzte2\nk5LQvYCdOxgcFsTd88aTW3aEVdu8ZyvNvYeP/WP7y6l9b3+p+kfDX1kmPDiARxZnUlLdxNMe1qt+\n8M0dHG5q48nrJxMSaF2753RLpiQxfXQMj+bsprbpuNXlDFhnl+G2VXkE+vvx6GKd3eNIGv7KUjNS\nus/mnttcSsHBBqvL6Zd3dlSy5suD/HDWODKHRVldzj8REVZcl8Gx4x08muP502n/96O95JYd4cFr\n00mMOrvtL9WZafgry917zQRiw4K4dVUebR3uvU5N3bE27lmTz4Qhkfxo1jiry+nRuPgIbr5sDKu3\nlfOJB0+nLa1p4vF3CpmdFs/iKUlWl+N1NPyV5aIGBbLiukx2Vzby3OY9VpdzRve9XkBDSztPfWMS\nQQHu+7/Pj2YlMzwmlHtfL3D7N9SedNpn94QE+vPIYp3d4wzu+69X+ZQrJySwYNJQfvN+MYWVjVaX\n06M38w7xdn4FP7sihbTESKvLOaPQIH+WLcigpLqJ328ptbqcs/bSh3vZVnaEBxdMICFS2z3OMKDw\nF5EYEXlXRIrtX3uc7yYinSLylf22diBjKu/14IJ0IkMCuX11Hh1utkxxdWMr971RwKTh0fzHZWOs\nLqdfZqXFMy8zkac3FLO/ttnqcvqtpLqJJ9YXcsX4BBZN1naPswz0zP9OYIMxJhnYYH/ckxZjzGT7\nbcEAx1ReKiYsiIcWppNX3sCLH+61upyTjDHc/fcCWto6efL6SQT4e84fzPfPTyfAT7h/bYFHbPre\naV+7JzTIn0cWZ2i7x4kG+q94IfCy/f7LwKIBvp7ycddkDuHq9ASefLeIPTVNVpcDwGtfHOS9XVXc\ndnUq4+LDrS7nrCRGhXDLValsKqxhXUGl1eX06YUtpXy5v56HFqQTH6HtHmcaaPgnGGMqAOxf43s5\nLkREckXkExHRNwjVKxHh4YUZhAb6c8fq7ZZfqVrR0MJDb+7g/FGDueni0ZbWcq6+c+FI0odG8uCb\nO9x6Ke2S6kaefLeIq9O7P/9RztVn+IvIeyJS0MNt4VmMM8IYkwV8C/iViIztZayb7W8SuTU13rnh\nh+pbfGQI98+fQG7ZEf64dZ9ldRhjuH31djo6DU9cPwl/D90cPMDfjxXXZVLdeJyn1hdZXU6POjq7\nuGXVdsKC/Fm+SGf3uEKf4W+MucIYk9HD7Q2gSkSGANi/VvfyGofsX0uBTcB5vRz3vDEmyxiTZbM5\nfyck5b4WT0liVqqNx9YVcqDOmg8rX/nsAFuKD3PXvDRGxoZZUoOjTB4ezbenj+QPH+91y4vpfr9l\nL3kH6lm2MANbRLDV5fiEgbZ91gLftd//LvDG6QeIyGARCbbfjwMuBnYOcFzl5USERxZnEuAn3PHa\ndpd/WHmgrpkVb+/korGxfHv6SJeO7Sy3Xp1KTFgw96zJp9ONFn4rqmrkl+8WMTcjkfkTh1hdjs8Y\naPivBK4UkWLgSvtjRCRLRF6wHzMeyBWRPGAjsNIYo+Gv+jQkKpS7rxnPx3tq+evnrluo7MRuUSLC\nL5ZOxM9D2z2niwoN5L7548krb+Avn+23uhygu91z66o8wkMCeHiRzu5xpQGFvzGm1hgz2xiTbP9a\nZ38+1xjz7/b7HxtjMo0xk+xfX3RE4co33HD+cC4eF8uKt3dxyEWblL+8dR+flNZx33zv2xx8waSh\nXDIujl+s2011Y6vV5fC7D0rZXt7AwwsziAvXdo8rec6EZeWTRISViyfS2WW4Z02+09s/pTVNJ5cP\n9sbdokSEhxdlcLyji+Vv7bK0lsLKRn71XhHXZA7hGm33uJyGv3J7w2MGccecVDYW1rDmy4NOG6ez\ny3DrqjyC/P1Y6cWbg4+OC+O/Z45lbd4hthRbM6uu3d7uiQwJZNnCdEtq8HUa/sojfOfCUWSNHMxD\nb+6k+qhz2hUvbCnli/31PLQw3evXk/nPGWMZHRfGfa8X0Nru+o10ntu0h/yDDSxflEGstnssoeGv\nPIKfX/eHr63tndz7uuOXKiiu6r7A6KoJvrGeTEigPw8vzGBfbTPPbnLtSqq7Ko7y9PvFXDtpKHMz\ntd1jFQ1/5THG2ML5+ZUprN9Zxdv5FQ573e4LjPIIC/JnxXW+c4HRJclxLJw8lGc37aHURUtpnGj3\nRIUG8tACbfdYScNfeZTvXzLdJb1NAAAJ/UlEQVSaScOieOCNHQ7bpvDZTXvYXt7Aiusyfe4Co3uu\nGU9woB/3veGahd/+Z+Medhw6yvJFmcSEBTl9PNU7DX/lUQL8/fjF0kkcbW3noTcHfrnIjkMN/HpD\ndwting+2IOIjQrh9ThofldSyNu+QU8facaiB37xfzMLJQ5mTkejUsVTfNPyVx0lNjODHlyezNu8Q\n63ec+0qVbR1d3PJqHoPDgljmwy2Ib00bwaTh0Tz81k4amtudMkZbRxe3rtpO9KAgHrzWd3/X7kTD\nX3mk/5o5lglDIrn39YJzDqynNxSzu7KRR6/LZLAPtyD8/YQVizKoO9bG4+uds+n7MxtL2FVxlEeu\ny/Dp37U70fBXHinQ349fLJ1I7bE2lr999u2fvAP1PLt5D0umDOOKCQlOqNCzZCRF8b2LRvPnT/fz\n5f4jDn3tgoMNPLOxhOvOS+KqdG33uAsNf+WxMpKi+K8ZY1m1rZzNRf2/WKm1vZNbVuVhCw/m/msn\nOLFCz/Lzq1JIiAjhnjUFDttGs7vdk0dMWBAP6O/arWj4K4/249njGBcfzl2vbaextX/tn6feLaKk\nuonHlk4kKjTQyRV6jvDgAB64dgI7K47y8tYyh7zmb963t9YWZxI9SNs97kTDX3m04AB/Hl86kcqj\nrTy2ru9+de6+On6/pZRvTR/BjBTdM+J0czISmZVq46n1hVQ0DGwhvfzyBv5nU3drbfZ4ba25Gw1/\n5fHOGzGY718ymj99sp+te2p7Pa65rYNbVuWRFB3K3fPGu7BCzyEiLFuYQUeXYdkAptIe7+jk1lV5\nxIUHaWvNTWn4K6/w8ytTGRU7iDte205zW8/71D6Ws5uy2mYeXzqJ8OAAF1foOYbHDOIns5PJKajk\n/d1V5/QaT28oprCqkZWLtbXmrjT8lVcIDfLnsSUT2V/XzJM97FP7cclhXt5axvcuGsWFY2MtqNCz\n/ODSMYyLD+f+N3bQ0nZ2C7/lHajnuc2lXD91GLPS4p1UoRooDX/lNaaPieU7F47kpY/2sq3sH9MV\nG1vbuW31dkbHhXHHnDQLK/QcQQF+rFiUQfmRFn7zfnG/f661vbvdYwsP5t752u5xZxr+yqvcPieN\noVGh3L467+RSxY9k76KioYUnrp9IaJC/xRV6juljYlk6dRjPf1BKUVVjv37m1xuKKa5uYuWSTG33\nuDkNf+VVwoMDWLkkkz01x/j1hmI2FVbzymcH+MGlY5g6Msbq8jzOXXPTCA8J4N41fS/89uX+I/xu\n8x6+kTWcmana7nF3Gv7K61yabOMbWcN5/oNSbnk1j+T4cP7flSlWl+WRYsODuWtuGp/tq2P1tvJe\njzvR7kmMDOGe+TqTyhNo+CuvdM/88djCg6lvaeepf5lMSKC2e87V9VOHkzVyMI9k7+LIsbYej/nl\ne0XsqTnGyiUTiQzRdo8n0PBXXikyJJD/+/40/nDT+WQOi7K6HI/m5ycsvy6DxtYOVuZ8/UK6L/Yf\n4fcflPLNacO5TC+c8xga/sprJSdEcGmyhpEjpCVG8v1LR/O33AN8vq/u5PMn2j1DovTCOU+j4a+U\n6pefzk4mKTqUe9bk025f+O3J9YWU1hzjsSUTidB2j0fR8FdK9cugoAAeWpBOUVUTL364l21ldbzw\n4V5unD6CS5LjrC5PnSW9xl0p1W9XTEjgqgkJ/Oq9IuIjQhgaFcpd2u7xSHrmr5Q6Kw8sSMdPhP11\nzTy+dKKuk+ShBhT+InK9iOwQkS4RyTrDcXNEpFBESkTkzoGMqZSyVlJ0KL/55nk8ujiTi8Zpu8dT\nDfQtuwBYDPyutwNExB94BrgSKAc+F5G1xphzXy9WKWUpXZ/f8w0o/I0xu6B7DfAzmAaUGGNK7cf+\nFVgIaPgrpZRFXNHzTwIOnPK43P7c14jIzSKSKyK5NTX935NVKaXU2enzzF9E3gMSe/jWPcaYN/ox\nRk9/FvS4QpQx5nngeYCsrKwzryKllFLqnPUZ/saYKwY4Rjkw/JTHw4BDA3xNpZRSA+CKts/nQLKI\njBaRIOAGYK0LxlVKKdWLgU71vE5EyoELgbdF5B3780NFJBvAGNMB/Ah4B9gFvGqM2TGwspVSSg3E\nQGf7rAHW9PD8IWDeKY+zgeyBjKWUUspx9ApfpZTyQdLX1mxWEZEaoGwALxEHHHZQOc7mSbWCZ9Xr\nSbWCZ9XrSbWCZ9U7kFpHGmP6XMvcbcN/oEQk1xjT65IT7sSTagXPqteTagXPqteTagXPqtcVtWrb\nRymlfJCGv1JK+SBvDv/nrS7gLHhSreBZ9XpSreBZ9XpSreBZ9Tq9Vq/t+SullOqdN5/5K6WU6oXX\nhb8nbRwjIi+JSLWIFFhdS19EZLiIbBSRXfYNfH5qdU1nIiIhIvKZiOTZ633I6pr6IiL+IvKliLxl\ndS19EZF9IpIvIl+JSK7V9ZyJiESLyGoR2W3/93uh1TX1RkRS7b/TE7ejIvIzp4zlTW0f+8YxRZyy\ncQzwTXfdOEZELgOagD8aYzKsrudMRGQIMMQY84WIRADbgEVu/LsVIMwY0yQigcCHwE+NMZ9YXFqv\nROTnQBYQaYyZb3U9ZyIi+4AsY4zbz5sXkZeBLcaYF+zriw0yxtRbXVdf7Hl2EJhujBnINU898rYz\n/5Mbxxhj2oATG8e4JWPMB0Cd1XX0hzGmwhjzhf1+I93rNPW4L4M7MN2a7A8D7Te3PdMRkWHANcAL\nVtfiTUQkErgMeBHAGNPmCcFvNxvY44zgB+8L/35vHKPOnYiMAs4DPrW2kjOzt1G+AqqBd40x7lzv\nr4DbgS6rC+knA6wXkW0icrPVxZzBGKAG+F97S+0FEQmzuqh+ugF4xVkv7m3h3++NY9S5EZFw4DXg\nZ8aYo1bXcybGmE5jzGS695CYJiJu2VoTkflAtTFmm9W1nIWLjTFTgLnAD+0tTHcUAEwBnjXGnAcc\nA9z6s0AAe3tqAbDKWWN4W/jrxjFOZO+dvwb82Rjzd6vr6S/7n/mbgDkWl9Kbi4EF9j76X4HLReRP\n1pZ0ZvaVezHGVNO9su80ayvqVTlQfspffavpfjNwd3OBL4wxVc4awNvCXzeOcRL7B6gvAruMMU9Z\nXU9fRMQmItH2+6HAFcBua6vqmTHmLmPMMGPMKLr/zb5vjPm2xWX1SkTC7B/6Y2+hXAW45Yw1Y0wl\ncEBEUu1PzQbccpLCab6JE1s+MMD1/N2NMaZDRE5sHOMPvOTOG8eIyCvATCDOvinOA8aYF62tqlcX\nA/8K5Nv76AB32/dqcEdDgJftMyb86N5EyO2nUHqIBGBN9/kAAcBfjDHrrC3pjH4M/Nl+QlgK3GRx\nPWckIoPonrH4H04dx5umeiqllOofb2v7KKWU6gcNf6WU8kEa/kop5YM0/JVSygdp+CullA/S8FdK\nKR+k4a+UUj5Iw18ppXzQ/wd0/h07Rm7kfgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(ws.oem_errors.value)# ?\n",
    "plt.plot(ws.y.value - ws.yf.value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([290.11153327,   5.92299149])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ws.x.value\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([290.2,   3.5])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ws.xa.value"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}