{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples for remorte sensing applications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "sys.path.insert(0, os.path.abspath(\"..\"))\n",
    "\n",
    "import geoclide as gc\n",
    "import math\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find the x and y components of the satellite position knowing its altitude and its viewing zenith and azimuth angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "vza = 45. # viewing zenith angle in degrees\n",
    "vaa = 45. # viewing azimuth angle in degrees\n",
    "sat_altitude = 700.  # satellite altitude in kilometers\n",
    "origin = gc.Point(0., 0., 0.) # origin is the viewer seeing the satellite\n",
    "# The vaa start from north going clockwise.\n",
    "# Let's assume that in our coordinate system the x axis is in the north direction\n",
    "# Then theta (zenith) angle = vza and phi (azimuth) angle = -vaa\n",
    "theta = vza\n",
    "phi = -vaa\n",
    "\n",
    "# Get the vector from ground to the satellite\n",
    "dir_to_sat = gc.ang2vec(theta=theta, phi=phi)\n",
    "ray = gc.Ray(o=origin, d=dir_to_sat) # create the ray, starting from origin going in dir_to_sat direction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here without considering the sphericity of the earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 494.97474683, -494.97474683,  700.        ])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b1 = gc.BBox(p1=gc.Point(-math.inf, -math.inf, 0.), p2=gc.Point(math.inf, math.inf, sat_altitude))\n",
    "ds_pp = gc.calc_intersection(b1, ray) # return an xarray dataset\n",
    "ds_pp['phit'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here with the consideration of the sphericity of the earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 472.61058011, -472.61058011,  668.37229212])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "earth_radius = 6378. # the equatorial earth radius in kilometers\n",
    "oTw = gc.get_translate_tf(gc.Vector(0., 0., -earth_radius))\n",
    "sphere_sat_alti = gc.Sphere(radius=earth_radius+sat_altitude, oTw=oTw)  # apply oTw to move the sphere center to earth center\n",
    "ds_sp = gc.calc_intersection(sphere_sat_alti, ray) # return an xarray dataset\n",
    "ds_sp['phit'].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Satellite camera directions (3MI example)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute all the pixel directions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of camera pixels in x and y\n",
    "nx = 509\n",
    "ny = 255\n",
    "\n",
    "# Satellite altitude (in km)\n",
    "z_sat = 830.\n",
    "\n",
    "# Nadir resolution (in km2)\n",
    "nad_res = 4.\n",
    "\n",
    "# swath (in km2)\n",
    "swath = 2200\n",
    "\n",
    "# Field of view of 1 pixel and of the camera\n",
    "pixel_fov = np.rad2deg(np.arctan(nad_res/z_sat))\n",
    "full_fov = np.rad2deg(np.arctan(0.5*swath/z_sat)) - pixel_fov\n",
    "\n",
    "focal_length = ((nx-1)*0.5) / ((np.tan(np.radians(full_fov))))\n",
    "focal_pos  = gc.Point(x=0., y=0., z=focal_length)\n",
    "\n",
    "x_ = -(nx-1)*0.5 + np.arange(nx)\n",
    "y_ = -(ny-1)*0.5 + np.arange(ny)\n",
    "x, y = np.meshgrid(x_, y_)\n",
    "x = x.flatten()\n",
    "y = y.flatten()\n",
    "z = np.zeros_like(x)\n",
    "id_pixels = gc.Point(x,y,z) # id = (0, 0, 0) corresponds to the pixel at the center of the camera\n",
    "dir_pixels = gc.normalize(id_pixels - focal_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.73904935, -0.36952468, -0.56325622],\n",
       "       [-0.73772467, -0.37032029, -0.56446895],\n",
       "       [-0.73639141, -0.37111789, -0.56568472],\n",
       "       ...,\n",
       "       [ 0.73639141, -0.37111789, -0.56568472],\n",
       "       [ 0.73772467, -0.37032029, -0.56446895],\n",
       "       [ 0.73904935, -0.36952468, -0.56325622]], shape=(509, 3))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# directions of first y pixels (x, y, z components)\n",
    "dir_pixels.to_numpy().reshape(ny,nx,3)[0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select only pixels that view a specific box zone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# box of size 50 km2 in x and y, and 10 km in z\n",
    "box = gc.BBox(p1=gc.Point(-25., -25., 0.), p2=gc.Point(25.,25.,10.))\n",
    "\n",
    "# satellite position/pixel positions. We duplicate to get same size as the number of pixels\n",
    "sat_pos = gc.Point(np.zeros_like(dir_pixels.x), np.zeros_like(dir_pixels.x),\n",
    "                   np.full_like(dir_pixels.x, z_sat))\n",
    "# create th rays\n",
    "r_sat = gc.Ray(sat_pos, dir_pixels)\n",
    "is_intersection = box.is_intersection(r_sat)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[-0.02581156, -0.02581156, -0.99933354],\n",
       "        [-0.02065172, -0.02581465, -0.99945341],\n",
       "        [-0.01549024, -0.02581706, -0.99954666],\n",
       "        [-0.01032751, -0.02581878, -0.99961329],\n",
       "        [-0.00516396, -0.02581982, -0.99965328],\n",
       "        [ 0.        , -0.02582016, -0.9996666 ],\n",
       "        [ 0.00516396, -0.02581982, -0.99965328],\n",
       "        [ 0.01032751, -0.02581878, -0.99961329],\n",
       "        [ 0.01549024, -0.02581706, -0.99954666],\n",
       "        [ 0.02065172, -0.02581465, -0.99945341]]),\n",
       " array([[-5., -5.,  0.],\n",
       "        [-4., -5.,  0.],\n",
       "        [-3., -5.,  0.],\n",
       "        [-2., -5.,  0.],\n",
       "        [-1., -5.,  0.],\n",
       "        [ 0., -5.,  0.],\n",
       "        [ 1., -5.,  0.],\n",
       "        [ 2., -5.,  0.],\n",
       "        [ 3., -5.,  0.],\n",
       "        [ 4., -5.,  0.]]))"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# pixels directions and x and y id of pixels\n",
    "dir_pixels.to_numpy()[is_intersection][0:10], id_pixels.to_numpy()[is_intersection][0:10]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "geoclide-py313",
   "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.13.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}