Proj4_Resample_UCERF_velo_05192021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# This code samples the UCERF3 velocity consensus estimates on the boundary \n",
    "Jey Kim (May 19, 2021)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as scipy\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "OUTPUT = open(\"vel_on_boundary.dat\",'w')\n",
    "#OUTPUT.write('%f %f %f %f %f %f %f\\n' % (lon,lat,vx,vy,sx,sy,cxy))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "243.32 243.46 33.86 33.96\n"
     ]
    }
   ],
   "source": [
    "#boundary lon, lat\n",
    "lon_min=np.round(-116.680+360,3)\n",
    "lon_max=np.round(-116.54+360,3)\n",
    "lat_min=33.860\n",
    "lat_max=33.960\n",
    "print(lon_min,lon_max,lat_min,lat_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt(\"vel_UCERF_NA.out\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon = data[:,0]\n",
    "lat = data[:,1]\n",
    "vx = data[:,2] # in [mm] \n",
    "vy = data[:,3] # The reference frame is NA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_sample = np.loadtxt(\"vel_001_100.gmt\")\n",
    "# boundary sampled at KNOTPOINTS (May 19, 2021)\n",
    "# JK"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon_bnd = boundary_sample[:,0]\n",
    "lat_bnd = boundary_sample[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11421 165\n"
     ]
    }
   ],
   "source": [
    "n_data = len(lon)\n",
    "n_bnd = len(lon_bnd)\n",
    "print(n_data,n_bnd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon_r=np.round(lon, 3)\n",
    "lat_r=np.round(lat, 3)\n",
    "lon_bnd_r=np.round(lon_bnd,3)\n",
    "lat_bnd_r=np.round(lat_bnd,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#for i in range(n_bnd):\n",
    "for i in range(n_bnd):\n",
    "#    print(lon_bnd_r[i],lat_bnd_r[i])\n",
    "    \n",
    "    if ((lon_bnd_r[i]== lon_min) or \\\n",
    "        (lon_bnd_r[i]== lon_max) or \\\n",
    "        (lat_bnd_r[i]== lat_min) or \\\n",
    "        (lat_bnd_r[i]== lat_max)):\n",
    "        \n",
    "        idx_target = np.argwhere((lon_r==lon_bnd_r[i]) & (lat_r==lat_bnd_r[i]))\n",
    "#         print(idx_target);\n",
    "#         print(lon_r[idx_target],lat_r[idx_target],vx[idx_target],vy[idx_target]);\n",
    "    \n",
    "    \n",
    "        lon_save = float(lon_r[idx_target])\n",
    "        lat_save = float(lat_r[idx_target])\n",
    "        vx_save = float(vx[idx_target])\n",
    "        vy_save = float(vy[idx_target])\n",
    "#         print(lon_save,lat_save,vx_save,vy_save);\n",
    "#         print();\n",
    "        OUTPUT.write('%f %f %f %f\\n' % (lon_save,lat_save,vx_save,vy_save))\n",
    "        \n",
    "    else: \n",
    "        continue\n",
    "#         print(\"pass inside\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#idx_lithosphere = np.argwhere((flag_pressure_np==1.0) | (flag_pressure_np==2.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}