Sunrise - Render Observation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "de49e17b",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Start\" data-toc-modified-id=\"Start-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Start</a></span></li><li><span><a href=\"#Setup\" data-toc-modified-id=\"Setup-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Setup</a></span><ul class=\"toc-item\"><li><span><a href=\"#Date-(year,-month,-day)-to-location-(lat,-lng)\" data-toc-modified-id=\"Date-(year,-month,-day)-to-location-(lat,-lng)-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Date (year, month, day) to location (lat, lng)</a></span></li><li><span><a href=\"#Location-(lat,-lng,-alt)-to-spatial-(x,-y,-z)\" data-toc-modified-id=\"Location-(lat,-lng,-alt)-to-spatial-(x,-y,-z)-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Location (lat, lng, alt) to spatial (x, y, z)</a></span></li><li><span><a href=\"#Tile-(zoom,-x,-y)-to-location-(lat,-lng)\" data-toc-modified-id=\"Tile-(zoom,-x,-y)-to-location-(lat,-lng)-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Tile (zoom, x, y) to location (lat, lng)</a></span></li></ul></li><li><span><a href=\"#Initialize-&amp;-Load-Data\" data-toc-modified-id=\"Initialize-&amp;-Load-Data-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Initialize &amp; Load Data</a></span></li><li><span><a href=\"#OSPRay\" data-toc-modified-id=\"OSPRay-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>OSPRay</a></span><ul class=\"toc-item\"><li><span><a href=\"#Config\" data-toc-modified-id=\"Config-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Config</a></span></li><li><span><a href=\"#Download\" data-toc-modified-id=\"Download-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>Download</a></span></li><li><span><a href=\"#Extract\" data-toc-modified-id=\"Extract-4.3\"><span class=\"toc-item-num\">4.3&nbsp;&nbsp;</span>Extract</a></span></li><li><span><a href=\"#Load\" data-toc-modified-id=\"Load-4.4\"><span class=\"toc-item-num\">4.4&nbsp;&nbsp;</span>Load</a></span></li><li><span><a href=\"#Init\" data-toc-modified-id=\"Init-4.5\"><span class=\"toc-item-num\">4.5&nbsp;&nbsp;</span>Init</a></span></li><li><span><a href=\"#Definitions\" data-toc-modified-id=\"Definitions-4.6\"><span class=\"toc-item-num\">4.6&nbsp;&nbsp;</span>Definitions</a></span></li></ul></li><li><span><a href=\"#Render\" data-toc-modified-id=\"Render-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Render</a></span><ul class=\"toc-item\"><li><span><a href=\"#Terrain-(OSPGeometry)\" data-toc-modified-id=\"Terrain-(OSPGeometry)-5.1\"><span class=\"toc-item-num\">5.1&nbsp;&nbsp;</span>Terrain (OSPGeometry)</a></span></li><li><span><a href=\"#Colormap-(OSPMaterial[])\" data-toc-modified-id=\"Colormap-(OSPMaterial[])-5.2\"><span class=\"toc-item-num\">5.2&nbsp;&nbsp;</span>Colormap (OSPMaterial[])</a></span></li><li><span><a href=\"#Observation-(uint[8])\" data-toc-modified-id=\"Observation-(uint[8])-5.3\"><span class=\"toc-item-num\">5.3&nbsp;&nbsp;</span>Observation (uint[8])</a></span></li><li><span><a href=\"#OSPWorld\" data-toc-modified-id=\"OSPWorld-5.4\"><span class=\"toc-item-num\">5.4&nbsp;&nbsp;</span>OSPWorld</a></span></li><li><span><a href=\"#OSPCamera\" data-toc-modified-id=\"OSPCamera-5.5\"><span class=\"toc-item-num\">5.5&nbsp;&nbsp;</span>OSPCamera</a></span></li><li><span><a href=\"#OSPRenderer\" data-toc-modified-id=\"OSPRenderer-5.6\"><span class=\"toc-item-num\">5.6&nbsp;&nbsp;</span>OSPRenderer</a></span></li><li><span><a href=\"#Load-Trail\" data-toc-modified-id=\"Load-Trail-5.7\"><span class=\"toc-item-num\">5.7&nbsp;&nbsp;</span>Load Trail</a></span></li><li><span><a href=\"#Render-Single\" data-toc-modified-id=\"Render-Single-5.8\"><span class=\"toc-item-num\">5.8&nbsp;&nbsp;</span>Render Single</a></span></li></ul></li><li><span><a href=\"#Video\" data-toc-modified-id=\"Video-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Video</a></span><ul class=\"toc-item\"><li><span><a href=\"#Animation-Widget\" data-toc-modified-id=\"Animation-Widget-6.1\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Animation Widget</a></span></li><li><span><a href=\"#ffmpeg\" data-toc-modified-id=\"ffmpeg-6.2\"><span class=\"toc-item-num\">6.2&nbsp;&nbsp;</span>ffmpeg</a></span></li><li><span><a href=\"#render()\" data-toc-modified-id=\"render()-6.3\"><span class=\"toc-item-num\">6.3&nbsp;&nbsp;</span><code>render()</code></a></span></li><li><span><a href=\"#Animate\" data-toc-modified-id=\"Animate-6.4\"><span class=\"toc-item-num\">6.4&nbsp;&nbsp;</span>Animate</a></span></li></ul></li><li><span><a href=\"#Figure-out-OSPRay-equirectangular-coordinates\" data-toc-modified-id=\"Figure-out-OSPRay-equirectangular-coordinates-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Figure out OSPRay equirectangular coordinates</a></span></li><li><span><a href=\"#Sunrise-as-a-Service\" data-toc-modified-id=\"Sunrise-as-a-Service-8\"><span class=\"toc-item-num\">8&nbsp;&nbsp;</span>Sunrise-as-a-Service</a></span><ul class=\"toc-item\"><li><span><a href=\"#Inputs\" data-toc-modified-id=\"Inputs-8.1\"><span class=\"toc-item-num\">8.1&nbsp;&nbsp;</span>Inputs</a></span></li><li><span><a href=\"#patchbay.pub\" data-toc-modified-id=\"patchbay.pub-8.2\"><span class=\"toc-item-num\">8.2&nbsp;&nbsp;</span><code>patchbay.pub</code></a></span></li><li><span><a href=\"#index.html\" data-toc-modified-id=\"index.html-8.3\"><span class=\"toc-item-num\">8.3&nbsp;&nbsp;</span><code>index.html</code></a></span></li><li><span><a href=\"#Run-Server\" data-toc-modified-id=\"Run-Server-8.4\"><span class=\"toc-item-num\">8.4&nbsp;&nbsp;</span>Run Server</a></span></li></ul></li><li><span><a href=\"#End\" data-toc-modified-id=\"End-9\"><span class=\"toc-item-num\">9&nbsp;&nbsp;</span>End</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2baae8f",
   "metadata": {},
   "source": [
    "# Start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99cfa5ab",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:45.141607Z",
     "start_time": "2023-08-29T13:05:44.890205Z"
    }
   },
   "outputs": [],
   "source": [
    "!hostname\n",
    "!date\n",
    "# %cd ~/src/Sunrise-Demo/\n",
    "from __future__ import annotations\n",
    "try:\n",
    "    from mediocreatbest import auto, run\n",
    "except ImportError:\n",
    "    %pip install --quiet --upgrade pip\n",
    "    %pip install --upgrade --force-reinstall \\\n",
    "        mediocreatbest@git+https://gist.github.com/player1537/3457b026ed6ef6696d758517f55a58df.git\n",
    "    from mediocreatbest import auto, run\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf47b5b5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:45.150643Z",
     "start_time": "2023-08-29T13:05:45.146005Z"
    }
   },
   "outputs": [],
   "source": [
    "DO_SAVE = True\n",
    "TERRAIN_QUALITY = (\n",
    "    0\n",
    "    # 1\n",
    "    # 2\n",
    ")\n",
    "VIDEO_QUALITY = (\n",
    "#     0\n",
    "#     1\n",
    "    2\n",
    ")\n",
    "\n",
    "# # Zoom level 13\n",
    "# RLO=824832\n",
    "# RHI=827647\n",
    "RLO=823296\n",
    "RHI=829438\n",
    "# CLO=559104\n",
    "CLO=558080\n",
    "# CHI=564479\n",
    "CHI=566270\n",
    "\n",
    "# # Zoom level 15\n",
    "RLO=3300096\n",
    "RHI=3310591\n",
    "CLO=2236928\n",
    "CHI=2259455\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73d798d3",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b863543",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "616ad2af",
   "metadata": {},
   "source": [
    "## Date (year, month, day) to location (lat, lng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fa1538e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:56.309017Z",
     "start_time": "2023-08-29T13:05:45.412670Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%scope -o convert_date_to_sun_location\n",
    "auto.skyfield; import skyfield.api\n",
    "\n",
    "planets = auto.skyfield.api.load('de421.bsp')\n",
    "sun = planets['sun']\n",
    "earth = planets['earth']\n",
    "\n",
    "ts = auto.skyfield.api.load.timescale()\n",
    "    \n",
    "def convert_date_to_sun_location(\n",
    "    now: auto.datetime.datetime,\n",
    "    /,\n",
    ") -> tuple[Degree, Degree]:\n",
    "    t = ts.from_datetime(now)\n",
    "\n",
    "    position = earth.at(t).observe(sun).apparent()\n",
    "\n",
    "    location = auto.skyfield.toposlib.wgs84.geographic_position_of(position)\n",
    "    latitude = location.latitude.degrees\n",
    "    longitude = location.longitude.degrees\n",
    "    \n",
    "    return latitude, longitude\n",
    "\n",
    "goal: 'seconds' = 10\n",
    "rate: '1/seconds' = 130.0\n",
    "N = int(goal * rate)\n",
    "\n",
    "keys = ['latitude', 'longitude']\n",
    "K = len(keys)\n",
    "\n",
    "X = auto.numpy.linspace(0.0, 90.0, num=N)\n",
    "Y = auto.numpy.zeros(\n",
    "    shape=(N, K),\n",
    "    dtype='f4',\n",
    ")\n",
    "\n",
    "ref = auto.datetime.datetime(year=2023, month=1, day=1, tzinfo=auto.datetime.timezone.utc)\n",
    "\n",
    "for n, x in enumerate(auto.tqdm.tqdm(X)):\n",
    "    values = convert_date_to_sun_location(ref + auto.datetime.timedelta(days=x))\n",
    "    for k in range(K):\n",
    "        Y[n, k] = values[k]\n",
    "\n",
    "%matplotlib inline\n",
    "fig = auto.matplotlib.pyplot.figure()\n",
    "\n",
    "for k, key in enumerate(keys):\n",
    "    ax = fig.add_subplot(K, 1, k+1)\n",
    "    ax.plot(X, Y[:, k], label=f'{key}')\n",
    "    ax.legend()\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4289844",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-15T21:06:47.748725Z",
     "start_time": "2023-08-15T21:06:47.744400Z"
    }
   },
   "source": [
    "## Location (lat, lng, alt) to spatial (x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbfc5438",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:56.317210Z",
     "start_time": "2023-08-29T13:05:56.311825Z"
    }
   },
   "outputs": [],
   "source": [
    "def convert_location_to_spatial(\n",
    "    *,\n",
    "    lat: Degree,\n",
    "    lng: Degree,\n",
    "    alt: Meter,\n",
    "    math=auto.math,\n",
    ") -> Tuple[Meter, Meter, Meter]:\n",
    "    # Thanks https://gis.stackexchange.com/a/4148\n",
    "    \n",
    "    #> Note that \"Lat/Lon/Alt\" is just another name for spherical coordinates, and \n",
    "    #> phi/theta/rho are just another name for latitude, longitude, and altitude.\n",
    "    #> :) (A minor difference: altitude is usually measured from the surface of the \n",
    "    #> sphere; rho is measured from the center -- to convert, just add/subtract the \n",
    "    #> radius of the sphere.)\n",
    "    phi: Radian = math.radians(lat)\n",
    "    theta: Radian = math.radians(lng)\n",
    "    \n",
    "    # Thanks https://en.wikipedia.org/wiki/Earth_radius\n",
    "    #> A globally-average value is usually considered to be 6,371 kilometres (3,959 mi)\n",
    "    rho: Meter = 6_371_100 + alt\n",
    "    # rho: Meter = 6_371_000 + alt\n",
    "    \n",
    "    #> x = math.cos(phi) * math.cos(theta) * rho\n",
    "    x: Meter = math.cos(phi) * math.cos(theta) * rho\n",
    "    \n",
    "    #> y = math.cos(phi) * math.sin(theta) * rho\n",
    "    y: Meter = math.cos(phi) * math.sin(theta) * rho\n",
    "\n",
    "    #> z = math.sin(phi) * rho # z is 'up'\n",
    "    z: Meter = math.sin(phi) * rho\n",
    "    \n",
    "    #> (Note there's some slightly arbitrary choices here in what each axis means...\n",
    "    #> you might want 'y' to point at the north pole instead of 'z', for example.)\n",
    "    \n",
    "    # I do :)\n",
    "    y, z = z, y\n",
    "    \n",
    "    return x, y, z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62591bd8",
   "metadata": {},
   "source": [
    "## Tile (zoom, x, y) to location (lat, lng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b1baf7e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:56.335813Z",
     "start_time": "2023-08-29T13:05:56.318971Z"
    }
   },
   "outputs": [],
   "source": [
    "# @auto.mediocreatbest.doctest\n",
    "def convert_tile_coordinates_to_latlng(\n",
    "    *,\n",
    "    z: int,\n",
    "    x: int,\n",
    "    y: int,\n",
    "    math=auto.numpy,\n",
    "    return_google_maps_url: bool=False,\n",
    ") -> Tuple[Degree, Degree]:\n",
    "    r\"\"\"\n",
    "\n",
    "    >>> convert_tile_coordinates_to_latlng(z=15, x=8739, y=12925, return_google_maps_url=True)\n",
    "    'https://www.google.com/maps/@35.48751102385376,-83.990478515625,15.79z?entry=ttu'\n",
    "\n",
    "    \"\"\"\n",
    "    \n",
    "    # Thanks https://gis.stackexchange.com/a/133535\n",
    "    \n",
    "    #> n = 2 ^ zoom\n",
    "    n = 2 ** z\n",
    "    \n",
    "    #> lon_deg = xtile / n * 360.0 - 180.0\n",
    "    lon_deg = x / n * 360.0 - 180.0\n",
    "    \n",
    "    #> lat_rad = arctan(sinh(π * (1 - 2 * ytile / n)))\n",
    "    lat_rad = math.arctan(math.sinh(math.pi * (1.0 - 2.0 * y / n)))\n",
    "    \n",
    "    #> lat_deg = lat_rad * 180.0 / π\n",
    "    lat_deg = lat_rad * 180.0 / math.pi\n",
    "    # lat_deg = -1.0 * lat_deg  # XXX(th): I think this is a bug in the original code\n",
    "    \n",
    "    if return_google_maps_url:\n",
    "        return f\"\"\"https://www.google.com/maps/@{lat_deg},{lon_deg},15.79z?entry=ttu\"\"\"\n",
    "    \n",
    "    return lat_deg, lon_deg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6ad621f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:05:56.345357Z",
     "start_time": "2023-08-29T13:05:56.338723Z"
    }
   },
   "outputs": [],
   "source": [
    "%%scope -o Tixel\n",
    "\n",
    "@auto.dataclasses.dataclass(eq=True, order=True, frozen=True)\n",
    "class Tixel:\n",
    "    y: int\n",
    "    x: int\n",
    "        \n",
    "    z: ClassVar[int] = 23\n",
    "        \n",
    "    @property\n",
    "    def location(self) -> Location:\n",
    "        \"\"\"Center lat/lng of current tixel\"\"\"\n",
    "        lat, lng = tile_coordinates_to_latlng(\n",
    "            z=1+self.z,\n",
    "            x=1+self.x,\n",
    "            y=1+self.y,\n",
    "        )\n",
    "        return Location(lat=lat, lng=lng)\n",
    "\n",
    "# @Tixel.register\n",
    "# def __tixel_from_coord(coord: Coord, /, *, res: int):\n",
    "#     dz = int(auto.math.log2(Tixel.z - coord.z))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c84508b",
   "metadata": {},
   "source": [
    "# Initialize & Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c03fcef",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:06:04.470357Z",
     "start_time": "2023-08-29T13:05:56.347079Z"
    }
   },
   "outputs": [],
   "source": [
    "@run\n",
    "def interpolate():\n",
    "    def load(\n",
    "        *,\n",
    "        name: str,\n",
    "        method: str,\n",
    "    ) -> auto.scipy.interpolate.RegularGridInterpolator:\n",
    "        path = auto.pathlib.Path.cwd() / 'temp' / f'park-{name}.npz'\n",
    "        assert path.exists(), \\\n",
    "            f\"\"\"{path=!r}\"\"\"\n",
    "        print(f\"\"\"Reading {path.stat().st_size :,d} bytes from {path}\"\"\")\n",
    "        \n",
    "        with auto.numpy.load(path) as f:\n",
    "            points0 = f['points0']\n",
    "            points1 = f['points1']\n",
    "            values = f['values']\n",
    "        \n",
    "        print(f\"\"\"  Interpolation ranges from:\n",
    "    {points0.shape} tixel rows\n",
    "    {points1.shape} tixel columns\"\"\")\n",
    "        \n",
    "        return auto.scipy.interpolate.RegularGridInterpolator(\n",
    "            points=(\n",
    "                points0,\n",
    "                points1,\n",
    "            ),\n",
    "            values=values,\n",
    "            bounds_error=False,\n",
    "            fill_value=0,\n",
    "            method=method,\n",
    "        )\n",
    "\n",
    "    elevation = load(\n",
    "        name='elevation',\n",
    "        method='linear',\n",
    "    )\n",
    "    \n",
    "    observation = load(\n",
    "        name='observation',\n",
    "        method='nearest',\n",
    "    )\n",
    "    \n",
    "    background = load(\n",
    "        name='background',\n",
    "        method='nearest',\n",
    "    )\n",
    "    \n",
    "    return auto.mediocreatbest.namedtuple.elevation.observation.background(\n",
    "        elevation,\n",
    "        observation,\n",
    "        background,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "461defcc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:07:54.330403Z",
     "start_time": "2023-08-29T13:07:49.348330Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='render', after=print)\n",
    "def terrain(interpolate, /):\n",
    "    \"\"\"mesh/OSPGeometry:mesh:position/vertex.position:vec3f[]+index/index:vec4ui[]\"\"\"\n",
    "    # rlo, rhi = 3300096, 3310591\n",
    "    # clo, chi = 2236928, 2259455\n",
    "\n",
    "    # real_rlo = interpolate.elevation.grid[0].min()\n",
    "    # real_rhi = interpolate.elevation.grid[0].max()\n",
    "    # real_clo = interpolate.elevation.grid[1].min()\n",
    "    # real_chi = interpolate.elevation.grid[1].max()\n",
    "    real_rlo = interpolate.observation.grid[0].min()\n",
    "    real_rhi = interpolate.observation.grid[0].max()\n",
    "    real_clo = interpolate.observation.grid[1].min()\n",
    "    real_chi = interpolate.observation.grid[1].max()\n",
    "    # real_rlo = interpolate.background.grid[0].min()\n",
    "    # real_rhi = interpolate.background.grid[0].max()\n",
    "    # real_clo = interpolate.background.grid[1].min()\n",
    "    # real_chi = interpolate.background.grid[1].max()\n",
    "    # assert rlo >= real_rlo, \\\n",
    "    #     f'{rlo=!r} < {real_rlo=!r}'\n",
    "    # assert rhi <= real_rhi, \\\n",
    "    #     f'{rhi=!r} > {real_rhi=!r}'\n",
    "    # assert clo >= real_clo, \\\n",
    "    #     f'{clo=!r} < {real_clo=!r}'\n",
    "    # assert chi <= real_chi, \\\n",
    "    #     f'{chi=!r} > {real_chi=!r}'\n",
    "\n",
    "    # rs = auto.numpy.arange(rlo, rhi)\n",
    "    # cs = auto.numpy.arange(clo, chi)\n",
    "    # rs = auto.numpy.arange(real_rlo, real_rhi)\n",
    "    # cs = auto.numpy.arange(real_clo, real_chi)\n",
    "    rs = auto.numpy.arange(RLO, RHI)\n",
    "    cs = auto.numpy.arange(CLO, CHI)\n",
    "    print(f'RS: {rs}')\n",
    "    print(f'CS: {cs}')\n",
    "\n",
    "    def decimate(\n",
    "        rs: auto.np.ndarray,\n",
    "        cs: auto.np.ndarray,\n",
    "        /,\n",
    "        *,\n",
    "        duration: 'seconds',\n",
    "        rate: 'item/s',\n",
    "    ) -> tuple[auto.np.ndarray, auto.np.ndarray]:\n",
    "        total_items = len(rs) * len(cs)\n",
    "\n",
    "        actual_items = duration * rate\n",
    "    #             print(f'{items = } = {duration} * {rate}')\n",
    "\n",
    "        factor_items = total_items / actual_items\n",
    "\n",
    "        stride = int(auto.math.ceil(auto.math.sqrt( factor_items )))\n",
    "\n",
    "    #             stride = int(auto.math.ceil())\n",
    "    #             print(f'{stride = } = int(auto.math.ceil( {len(rs)} / {items} == {len(rs)/items} ))')\n",
    "\n",
    "        rs = rs[::stride]\n",
    "        cs = cs[::stride]\n",
    "\n",
    "        print(f'Decimate by factor of {stride}**2={stride**2 :,d} from {total_items :,d} to {len(rs)*len(cs) :,d} vertices')\n",
    "        return rs, cs\n",
    "\n",
    "    rs, cs = decimate(\n",
    "        rs, cs,\n",
    "        duration=[\n",
    "            1,  # lo resolution\n",
    "            10,  # medium resolution\n",
    "            1000,  # high resolution\n",
    "        ][TERRAIN_QUALITY],\n",
    "        rate=6260094,\n",
    "    )\n",
    "\n",
    "    position: 'OSPGeometry:mesh:vertex.position:vec3f[]' = auto.numpy.zeros(\n",
    "        shape=(len(rs), len(cs)),\n",
    "        dtype=[\n",
    "            ('x', 'f4'),\n",
    "            ('y', 'f4'),\n",
    "            ('z', 'f4'),\n",
    "        ],\n",
    "    )\n",
    "    \n",
    "    texcoord: 'OSPGeometry:mesh:vertex.texcoord:vec2f[]' = auto.numpy.zeros(\n",
    "        shape=(len(rs), len(cs)),\n",
    "        dtype=[\n",
    "            ('u', 'f4'),\n",
    "            ('v', 'f4'),\n",
    "        ],\n",
    "    )\n",
    "        \n",
    "    normal: 'OSPGeometry:mesh:vertex.normal:vec3f[]' = auto.numpy.zeros(\n",
    "        shape=(len(rs), len(cs)),\n",
    "        dtype=[ ('x', 'f4'), ('y', 'f4'), ('z', 'f4') ],\n",
    "    )\n",
    "    \n",
    "    @auto.mediocreatbest.immediate\n",
    "    def __vert():\n",
    "        with auto.tqdm.auto.tqdm(total=position.size, unit='vert') as tqdm:\n",
    "            for ri, r in enumerate(rs):\n",
    "                lat, lngs = convert_tile_coordinates_to_latlng(\n",
    "                    z=Tixel.z,\n",
    "                    y=r,\n",
    "                    x=cs,\n",
    "                    math=auto.numpy,\n",
    "                )\n",
    "                # lat *= -1\n",
    "                lngs *= -1\n",
    "\n",
    "                if ri == 0:\n",
    "                    print(lat, lngs[0])\n",
    "                elif ri == len(rs) - 1:\n",
    "                    print(lat, lngs[-1])\n",
    "                \n",
    "                # r: float\n",
    "                # cs: list[float]\n",
    "                \n",
    "                # rs: list[float]\n",
    "                \n",
    "                # grid: list[tuple[float, float]]\n",
    "\n",
    "                grid = auto.numpy.hstack(tuple(\n",
    "                    v.reshape((-1, 1))\n",
    "                    for v in auto.numpy.meshgrid(\n",
    "                        [r],\n",
    "                        cs,\n",
    "                        indexing='ij',\n",
    "                    )\n",
    "                ))\n",
    "                assert grid.shape == (len(cs), 2), \\\n",
    "                    f\"\"\"Wrong shape: {grid.shape=!r}\"\"\"\n",
    "                # interpolate.elevation.fill_value = 0\n",
    "                alts = interpolate.elevation(\n",
    "                    grid,\n",
    "                    method='linear',\n",
    "                )\n",
    "                assert alts.shape == (len(cs),), \\\n",
    "                    f\"\"\"Wrong shape: {alts.shape=!r}\"\"\"\n",
    "\n",
    "                Zmin = 0; Zmax = 1\n",
    "                # # Zmin = 0.2527998\n",
    "                # # Zmax = 2.02335\n",
    "\n",
    "                # # xs = auto.numpy.linspace(0.0, 1.0, num=len(cs))\n",
    "                # # ys = auto.numpy.array([ri / len(rs)] * len(cs))\n",
    "                # xs = cs / (2 ** Tixel.z)\n",
    "                # xs = 1 - xs  # flip x\n",
    "                # # print(f'{cs[0]=!r} {Tixel.z=!r} {(Tixel.z**2)=!r}')\n",
    "                # assert auto.np.all(0.0 <= xs)\n",
    "                # assert auto.np.all(xs <= 1.0)\n",
    "                # ys = auto.numpy.array([r / (2 ** Tixel.z)] * len(cs))\n",
    "                # ys = 1 - ys  # flip y\n",
    "                # # print(f'{ri=}, {len(rs)=}, {Tixel.z=!r} {(Tixel.z**2)=!r}')\n",
    "                # assert auto.np.all(0.0 <= ys)\n",
    "                # assert auto.np.all(ys <= 1.0)\n",
    "                # zs = alts / 1000.0\n",
    "                # zs = (zs - Zmin) / (Zmax - Zmin)\n",
    "                # zs = -zs  # flip z\n",
    "                # assert xs.shape == (len(cs),), \\\n",
    "                #     f\"\"\"Wrong shape: {xs.shape=!r} {len(cs)=!r}\"\"\"\n",
    "                # assert ys.shape == (len(cs),), \\\n",
    "                #     f\"\"\"Wrong shape: {ys.shape=!r} {len(cs)=!r}\"\"\"\n",
    "                # assert zs.shape == (len(cs),), \\\n",
    "                #     f\"\"\"Wrong shape: {zs.shape=!r} {len(cs)=!r}\"\"\"\n",
    "\n",
    "                xs, ys, zs = convert_location_to_spatial(\n",
    "                    lat=lat,\n",
    "                    lng=lngs,\n",
    "                    alt=alts,\n",
    "                    math=auto.numpy,\n",
    "                )\n",
    "                # xs, ys, zs = (v/1000.0 for v in [xs, ys, zs])\n",
    "                # assert xs.shape == (len(cs),), \\\n",
    "                #     f\"\"\"Wrong shape: {xs.shape=!r}\"\"\"\n",
    "\n",
    "                position[ri, :]['x'] = xs\n",
    "                position[ri, :]['y'] = ys\n",
    "                position[ri, :]['z'] = zs\n",
    "                \n",
    "                normal[ri, :]['x'] = xs\n",
    "                normal[ri, :]['y'] = ys\n",
    "                normal[ri, :]['z'] = zs\n",
    "                \n",
    "                us = auto.numpy.linspace(0.0, 1.0, num=len(cs))\n",
    "                vs = ri / len(rs)\n",
    "                \n",
    "                texcoord[ri, :]['u'] = us\n",
    "                texcoord[ri, :]['v'] = vs\n",
    "\n",
    "                tqdm.update(len(cs))\n",
    "\n",
    "    print(f'Xmin={position[\"x\"][0,0]}, Xmax={position[\"x\"][-1,-1]}')\n",
    "    print(f'Ymin={position[\"y\"][0,0]}, Ymax={position[\"y\"][-1,-1]}')\n",
    "    print(f'Zmin={position[\"z\"][0,0]}, Zmax={position[\"z\"][-1,-1]}')\n",
    "\n",
    "    print(f'{position[\"x\"].min()=!r} ;; {position[\"x\"].max()=!r}')\n",
    "    print(f'{position[\"y\"].min()=!r} ;; {position[\"y\"].max()=!r}')\n",
    "    print(f'{position[\"z\"].min()=!r} ;; {position[\"z\"].max()=!r}')\n",
    "\n",
    "    print(f'{texcoord[\"u\"].min()=!r} ;; {texcoord[\"u\"].max()=!r}')\n",
    "    print(f'{texcoord[\"v\"].min()=!r} ;; {texcoord[\"v\"].max()=!r}')\n",
    "\n",
    "    # index = OSPGeometry:mesh:index:vec4ui[]\n",
    "    index = auto.numpy.zeros(\n",
    "        shape=(len(rs)-1, len(cs)-1),\n",
    "        dtype=[\n",
    "            ('a', 'u4'),\n",
    "            ('b', 'u4'),\n",
    "            ('c', 'u4'),\n",
    "            ('d', 'u4'),\n",
    "        ],\n",
    "    )\n",
    "\n",
    "    @auto.mediocreatbest.immediate\n",
    "    def __face():\n",
    "        def indexof(ri: int, ci: int, /) -> int:\n",
    "            return ri * len(cs) + ci\n",
    "\n",
    "        with auto.tqdm.auto.tqdm(total=index.size, unit='face') as tqdm:\n",
    "            # a -1- b\n",
    "            # |     |\n",
    "            # 4     2\n",
    "            # |     |\n",
    "            # d -3- c\n",
    "\n",
    "            for ri in range(0, len(rs)-1):\n",
    "                cis = auto.numpy.arange(len(cs)-1, dtype='u4')\n",
    "\n",
    "                a = indexof(ri+0, cis+0)\n",
    "                b = indexof(ri+0, cis+1)\n",
    "                c = indexof(ri+1, cis+1)\n",
    "                d = indexof(ri+1, cis+0)\n",
    "\n",
    "                index[ri, :]['a'] = a\n",
    "                index[ri, :]['b'] = b\n",
    "                index[ri, :]['c'] = c\n",
    "                index[ri, :]['d'] = d\n",
    "                tqdm.update(len(cis))\n",
    "\n",
    "        return index\n",
    "\n",
    "    print(f'{index[\"a\"].min()=!r} ;; {index[\"a\"].max()=!r}')\n",
    "    print(f'{index[\"b\"].min()=!r} ;; {index[\"b\"].max()=!r}')\n",
    "    print(f'{index[\"c\"].min()=!r} ;; {index[\"c\"].max()=!r}')\n",
    "    print(f'{index[\"d\"].min()=!r} ;; {index[\"d\"].max()=!r}')\n",
    "    \n",
    "    \n",
    "    def Flatten(data: auto.numpy.NDArray, /) -> auto.numpy.NDArray:\n",
    "        return data.reshape((-1,))\n",
    "    \n",
    "    position = Flatten(position)\n",
    "    normal = Flatten(normal)\n",
    "    texcoord = Flatten(texcoord)\n",
    "    index = Flatten(index)\n",
    "    \n",
    "    \n",
    "    #--- Save to Disk\n",
    "    \n",
    "    def Save(data: auto.numpy.NDArray, /, *, path: auto.pathlib.Path):\n",
    "        if not DO_SAVE:\n",
    "            return\n",
    "        \n",
    "        path.parent.mkdir(parents=True, exist_ok=True)\n",
    "        with open(path, 'wb') as f:\n",
    "            print(f'Writing {len(data.shape)}I = {data.shape=!r} to {path=!r}')\n",
    "            f.write(auto.struct.pack('I', len(data.shape)))\n",
    "            f.write(auto.struct.pack(f'{len(data.shape)}I', *data.shape))\n",
    "            data.tofile(f)\n",
    "        print(f\"Wrote {path.stat().st_size :,d} bytes to {path}\")\n",
    "    \n",
    "    Save(\n",
    "        position,\n",
    "        path=auto.pathlib.Path.cwd().joinpath(\n",
    "            f\"data\",\n",
    "            f\"park\",\n",
    "            f\"OSPGeometry.mesh.vertex.position.vec3f.bin\",\n",
    "        ),\n",
    "    )\n",
    "    \n",
    "    Save(\n",
    "        normal,\n",
    "        path=auto.pathlib.Path.cwd().joinpath(\n",
    "            f\"data\",\n",
    "            f\"park\",\n",
    "            f\"OSPGeometry.mesh.vertex.normal.vec3f.bin\",\n",
    "        ),\n",
    "    )\n",
    "    \n",
    "    Save(\n",
    "        texcoord,\n",
    "        path=auto.pathlib.Path.cwd().joinpath(\n",
    "            f\"data\",\n",
    "            f\"park\",\n",
    "            f\"OSPGeometry.mesh.vertex.texcoord.vec2f.bin\",\n",
    "        ),\n",
    "    )\n",
    "    \n",
    "    Save(\n",
    "        index,\n",
    "        path=auto.pathlib.Path.cwd().joinpath(\n",
    "            f\"data\",\n",
    "            f\"park\",\n",
    "            f\"OSPGeometry.mesh.index.vec4ui.bin\",\n",
    "        ),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "617df776",
   "metadata": {},
   "source": [
    "## Colormap (OSPMaterial[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d489d157",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:07:12.287682Z",
     "start_time": "2023-08-29T13:06:09.253945Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='render', after=print)\n",
    "def colormap(interpolate, /):\n",
    "    def Data(array: auto.numpy.ndarray, /, *, type: lib.OSPDataType) -> lib.OSPData:\n",
    "        if len(array.shape) == 0:\n",
    "            array = array[None, None, None]\n",
    "\n",
    "        elif len(array.shape) == 1:\n",
    "            array = array[:, None, None]\n",
    "\n",
    "        elif len(array.shape) == 2:\n",
    "            array = array[:, :, None]\n",
    "            \n",
    "        elif len(array.shape) == 3:\n",
    "            array = array[:, :, :]\n",
    "\n",
    "        else:\n",
    "            raise NotImplementedError()\n",
    "\n",
    "        src = lib.ospNewSharedData(\n",
    "            array.ctypes.data, type,\n",
    "            array.shape[0], array.strides[0],\n",
    "            array.shape[1], array.strides[1],\n",
    "            array.shape[2], array.strides[2],\n",
    "        )\n",
    "        lib.ospCommit(src)\n",
    "\n",
    "        dst = lib.ospNewData(type, *array.shape)\n",
    "        lib.ospCopyData(src, dst, 0, 0, 0)\n",
    "        lib.ospCommit(dst)\n",
    "\n",
    "        lib.ospRelease(src)\n",
    "        return dst\n",
    "\n",
    "    # Thanks https://species.atlasofthesmokies.org/style.css\n",
    "    pinks = auto.numpy.array([\n",
    "        # Black, obviously.\n",
    "        # Represents 0-25%\n",
    "        [0, 0, 0],\n",
    "\n",
    "        #> #legend-pink-50 { background-color: #f069a3; }\n",
    "        # Represents 25-50%\n",
    "        [247, 105, 164],\n",
    "\n",
    "        #> #legend-pink-75 { background-color: #c31d8e; }\n",
    "        # Represents 50-75%\n",
    "        [196, 28, 142],\n",
    "\n",
    "        #> #legend-pink-100 { background-color: #752879; }\n",
    "        # Represents 75-100%\n",
    "        [121, 0, 124],\n",
    "    ], dtype=auto.numpy.uint8)\n",
    "\n",
    "    # rlo, rhi = 3300096, 3310591\n",
    "    # clo, chi = 2236928, 2259455\n",
    "    # rlo = interpolate.background.grid[0].min()\n",
    "    # rhi = interpolate.background.grid[0].max()\n",
    "    # clo = interpolate.background.grid[1].min()\n",
    "    # chi = interpolate.background.grid[1].max()\n",
    "    # rlo = interpolate.elevation.grid[0].min()\n",
    "    # rhi = interpolate.elevation.grid[0].max()\n",
    "    # clo = interpolate.elevation.grid[1].min()\n",
    "    # chi = interpolate.elevation.grid[1].max()\n",
    "    rlo = interpolate.observation.grid[0].min()\n",
    "    rhi = interpolate.observation.grid[0].max()\n",
    "    clo = interpolate.observation.grid[1].min()\n",
    "    chi = interpolate.observation.grid[1].max()\n",
    "\n",
    "    print(f'elevation RLO: {interpolate.elevation.grid[0].min()}')\n",
    "    print(f'elevation CHI: {interpolate.elevation.grid[0].max()}')\n",
    "    print(f'elevation CLO: {interpolate.elevation.grid[1].min()}')\n",
    "    print(f'elevation CHI: {interpolate.elevation.grid[1].max()}')\n",
    "    print(f'background RLO: {interpolate.background.grid[0].min()}')\n",
    "    print(f'background CHI: {interpolate.background.grid[0].max()}')\n",
    "    print(f'background CLO: {interpolate.background.grid[1].min()}')\n",
    "    print(f'background CHI: {interpolate.background.grid[1].max()}')\n",
    "    print(f'observation RLO: {interpolate.observation.grid[0].min()}')\n",
    "    print(f'observation CHI: {interpolate.observation.grid[0].max()}')\n",
    "    print(f'observation CLO: {interpolate.observation.grid[1].min()}')\n",
    "    print(f'observation CHI: {interpolate.observation.grid[1].max()}')\n",
    "\n",
    "    rlo = RLO\n",
    "    rhi = RHI\n",
    "    clo = CLO\n",
    "    chi = CHI\n",
    "    # assert rlo >= real_rlo, \\\n",
    "    #     f'{rlo=!r} < {real_rlo=!r}'\n",
    "    # assert rhi <= real_rhi, \\\n",
    "    #     f'{rhi=!r} > {real_rhi=!r}'\n",
    "    # assert clo >= real_clo, \\\n",
    "    #     f'{clo=!r} < {real_clo=!r}'\n",
    "    # assert chi <= real_chi, \\\n",
    "    #     f'{chi=!r} > {real_chi=!r}'\n",
    "    rs = auto.numpy.linspace(rlo, rhi, num=1024 * 4)\n",
    "    rs = rs[::-1]  # Flip Y\n",
    "    cs = auto.numpy.linspace(clo, chi, num=2048 * 4)\n",
    "    # rs = auto.numpy.linspace(RLO, RHI, num=1024 * 4)\n",
    "    # rs = rs[::-1]  # Flip Y\n",
    "    # cs = auto.numpy.linspace(CLO, CHI, num=2048 * 4)\n",
    "    \n",
    "    @auto.mediocreatbest.immediate\n",
    "    @auto.mediocreatbest.coroutine\n",
    "    def display():\n",
    "        handle = auto.IPython.display.display(None, display_id=True)\n",
    "        \n",
    "        while True:\n",
    "            obj = yield\n",
    "            handle.update(obj)\n",
    "    \n",
    "    def Colormap(pi: int, /) -> lib.OSPTexture:\n",
    "        rgb = auto.numpy.zeros(\n",
    "            shape=(len(rs), len(cs)),\n",
    "            dtype=[\n",
    "                ('r', 'f4'),\n",
    "                ('g', 'f4'),\n",
    "                ('b', 'f4'),\n",
    "            ],\n",
    "        )\n",
    "        \n",
    "        for ri, r in enumerate(rs):\n",
    "            grid = auto.numpy.hstack(tuple(\n",
    "                v.reshape((-1, 1))\n",
    "                for v in auto.numpy.meshgrid(\n",
    "                    [r],\n",
    "                    cs,\n",
    "                    indexing='ij',\n",
    "                )\n",
    "            ))\n",
    "            assert grid.shape == (len(cs), 2), \\\n",
    "                f\"\"\"Wrong shape: {grid.shape=!r}\"\"\"\n",
    "#             print(f\"\"\"{grid[:, 0].min()=!r} ;; {grid[:, 0].max()=!r}\n",
    "# {grid[:, 1].min()=!r} ;; {grid[:, 1].max()=!r}\"\"\")\n",
    "            \n",
    "            rgbs = interpolate.background(\n",
    "                grid,\n",
    "                method='linear',\n",
    "            )\n",
    "            assert rgbs.shape == (len(cs), 3), \\\n",
    "                f\"\"\"Wrong shape: {rgbs.shape=!r}\"\"\"\n",
    "            \n",
    "            def composite(background: 'RGB', foreground: 'RGB', *, alpha: float):\n",
    "                return background * (1.0 - alpha) + foreground * alpha\n",
    "            \n",
    "            pink = pinks[None, pi, :]\n",
    "            pink = pink.astype(auto.numpy.float32)\n",
    "            pink /= 255.0\n",
    "            \n",
    "            if pi == 0:\n",
    "                alpha = 0.0\n",
    "            else:\n",
    "                alpha = 0.5\n",
    "            \n",
    "            rgbs = composite(rgbs, pink, alpha=alpha)\n",
    "            assert rgbs.shape == (len(cs), 3), \\\n",
    "                f\"\"\"Wrong shape: {rgbs.shape=!r}\"\"\"\n",
    "            \n",
    "            rgb['r'][ri, :] = rgbs[:, 0]\n",
    "            rgb['g'][ri, :] = rgbs[:, 1]\n",
    "            rgb['b'][ri, :] = rgbs[:, 2]\n",
    "\n",
    "#         import numpy.lib.recfunctions\n",
    "#         display.send(auto.PIL.Image.fromarray(\n",
    "#             (255.0 * auto.numpy.lib.recfunctions.structured_to_unstructured(rgb)).astype(auto.numpy.uint8),\n",
    "#             mode='RGB',\n",
    "#         ))\n",
    "\n",
    "\n",
    "        #--- Save to Disk\n",
    "\n",
    "        def Save(data: auto.numpy.NDArray, /, *, path: auto.pathlib.Path):\n",
    "            if not DO_SAVE:\n",
    "                return\n",
    "            \n",
    "            path.parent.mkdir(parents=True, exist_ok=True)\n",
    "            with open(path, 'wb') as f:\n",
    "                f.write(auto.struct.pack('I', len(data.shape)))\n",
    "                f.write(auto.struct.pack(f'{len(data.shape)}I', *data.shape))\n",
    "                data.tofile(f)\n",
    "            print(f\"Wrote {path.stat().st_size :,d} bytes to {path}\")\n",
    "\n",
    "        Save(\n",
    "            rgb,\n",
    "            path=auto.pathlib.Path.cwd().joinpath(\n",
    "                f\"data\",\n",
    "                f\"pink{pi}\",\n",
    "                f\"OSPTexture.texture2d.data.vec3f.bin\",\n",
    "            ),\n",
    "        )\n",
    "    \n",
    "    for pi in range(len(pinks)):\n",
    "        Colormap(pi)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ea5c450",
   "metadata": {},
   "source": [
    "## Observation (uint[8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57ab50eb",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T13:07:14.010892Z",
     "start_time": "2023-08-29T13:07:12.290858Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='render', after=print)\n",
    "def observation(interpolate, /) -> 'uint8[]':\n",
    "    \"\"\"OSPGeometricModel::index:uint8[]\"\"\"\n",
    "    \n",
    "    def Data(array: auto.numpy.ndarray, /, *, type: lib.OSPDataType) -> lib.OSPData:\n",
    "        if len(array.shape) == 0:\n",
    "            array = array[None, None, None]\n",
    "\n",
    "        elif len(array.shape) == 1:\n",
    "            array = array[:, None, None]\n",
    "\n",
    "        elif len(array.shape) == 2:\n",
    "            array = array[:, :, None]\n",
    "            \n",
    "        elif len(array.shape) == 3:\n",
    "            array = array[:, :, :]\n",
    "\n",
    "        else:\n",
    "            raise NotImplementedError()\n",
    "\n",
    "        src = lib.ospNewSharedData(\n",
    "            array.ctypes.data, type,\n",
    "            array.shape[0], array.strides[0],\n",
    "            array.shape[1], array.strides[1],\n",
    "            array.shape[2], array.strides[2],\n",
    "        )\n",
    "        lib.ospCommit(src)\n",
    "\n",
    "        dst = lib.ospNewData(type, *array.shape)\n",
    "        lib.ospCopyData(src, dst, 0, 0, 0)\n",
    "        lib.ospCommit(dst)\n",
    "\n",
    "        lib.ospRelease(src)\n",
    "        return dst\n",
    "    \n",
    "    # rlo, rhi = 3300096, 3310591\n",
    "    # clo, chi = 2236928, 2259455\n",
    "\n",
    "    real_rlo = interpolate.observation.grid[0].min()\n",
    "    real_rhi = interpolate.observation.grid[0].max()\n",
    "    real_clo = interpolate.observation.grid[1].min()\n",
    "    real_chi = interpolate.observation.grid[1].max()\n",
    "    # real_rlo = interpolate.background.grid[0].min()\n",
    "    # real_rhi = interpolate.background.grid[0].max()\n",
    "    # real_clo = interpolate.background.grid[1].min()\n",
    "    # real_chi = interpolate.background.grid[1].max()\n",
    "    # # real_rlo = interpolate.elevation.grid[0].min()\n",
    "    # # real_rhi = interpolate.elevation.grid[0].max()\n",
    "    # # real_clo = interpolate.elevation.grid[1].min()\n",
    "    # # real_chi = interpolate.elevation.grid[1].max()\n",
    "    # assert rlo >= real_rlo, \\\n",
    "    #     f'{rlo=!r} < {real_rlo=!r}'\n",
    "    # assert rhi <= real_rhi, \\\n",
    "    #     f'{rhi=!r} > {real_rhi=!r}'\n",
    "    # assert clo >= real_clo, \\\n",
    "    #     f'{clo=!r} < {real_clo=!r}'\n",
    "    # assert chi <= real_chi, \\\n",
    "    #     f'{chi=!r} > {real_chi=!r}'\n",
    "\n",
    "    rs = auto.numpy.arange(RLO, RHI)\n",
    "    cs = auto.numpy.arange(CLO, CHI)\n",
    "    # rs = auto.numpy.arange(real_rlo, real_rhi)\n",
    "    # cs = auto.numpy.arange(real_clo, real_chi)\n",
    "    # rs = auto.numpy.arange(rlo, rhi)\n",
    "    # cs = auto.numpy.arange(clo, chi)\n",
    "\n",
    "    def decimate(\n",
    "        rs: auto.np.ndarray,\n",
    "        cs: auto.np.ndarray,\n",
    "        /,\n",
    "        *,\n",
    "        duration: 'seconds',\n",
    "        rate: 'item/s',\n",
    "    ) -> tuple[auto.np.ndarray, auto.np.ndarray]:\n",
    "        total_items = len(rs) * len(cs)\n",
    "\n",
    "        actual_items = duration * rate\n",
    "    #             print(f'{items = } = {duration} * {rate}')\n",
    "\n",
    "        factor_items = total_items / actual_items\n",
    "\n",
    "        stride = int(auto.math.ceil(auto.math.sqrt( factor_items )))\n",
    "\n",
    "    #             stride = int(auto.math.ceil())\n",
    "    #             print(f'{stride = } = int(auto.math.ceil( {len(rs)} / {items} == {len(rs)/items} ))')\n",
    "\n",
    "        rs = rs[::stride]\n",
    "        cs = cs[::stride]\n",
    "\n",
    "        print(f'Decimate by factor of {stride}**2={stride**2 :,d} from {total_items :,d} to {len(rs)*len(cs) :,d} vertices')\n",
    "        return rs, cs\n",
    "\n",
    "    rs, cs = decimate(\n",
    "        rs, cs,\n",
    "        duration=[\n",
    "            1,  # lo resolution\n",
    "            10,  # medium resolution\n",
    "            1000,  # high resolution\n",
    "        ][TERRAIN_QUALITY],\n",
    "        rate=6260094,\n",
    "    )\n",
    "    \n",
    "    index = auto.numpy.zeros(\n",
    "        shape=(len(rs), len(cs)),\n",
    "        dtype='u1',\n",
    "    )\n",
    "    \n",
    "    for (ri, r) in enumerate(auto.tqdm.autonotebook.tqdm(rs, unit='row')):\n",
    "        grid = auto.numpy.hstack(tuple(\n",
    "            v.reshape((-1, 1))\n",
    "            for v in auto.numpy.meshgrid(\n",
    "                [r],\n",
    "                cs,\n",
    "                indexing='ij',\n",
    "            )\n",
    "        ))\n",
    "        assert grid.shape == (len(cs), 2), \\\n",
    "            f\"\"\"Wrong shape: {grid.shape=!r}\"\"\"\n",
    "#             print(f\"\"\"{grid[:, 0].min()=!r} ;; {grid[:, 0].max()=!r}\n",
    "# {grid[:, 1].min()=!r} ;; {grid[:, 1].max()=!r}\"\"\")\n",
    "        \n",
    "        try:\n",
    "            interpolate.observation.fill_value = 0\n",
    "            category = interpolate.observation(\n",
    "                grid,\n",
    "                method='nearest',\n",
    "            )\n",
    "        except ValueError:\n",
    "            print(f'{grid.shape=!r}')\n",
    "            print(f'{auto.np.min(grid, axis=0)=!r}')\n",
    "            print(f'{auto.np.max(grid, axis=0)=!r}')\n",
    "            raise\n",
    "\n",
    "        assert category.shape == (len(cs),), \\\n",
    "            f\"\"\"Wrong shape: {category.shape=!r}\"\"\"\n",
    "        assert auto.numpy.all((0 <= category) & (category < 4))\n",
    "        \n",
    "        index[ri, :] = category\n",
    "    \n",
    "    index = index.reshape((-1,))\n",
    "    print(auto.collections.Counter(index.flat))\n",
    "    \n",
    "    \n",
    "    #--- Save to Disk\n",
    "\n",
    "    def Save(data: auto.numpy.NDArray, /, *, path: auto.pathlib.Path):\n",
    "        if not DO_SAVE:\n",
    "            return\n",
    "        \n",
    "        path.parent.mkdir(parents=True, exist_ok=True)\n",
    "        with open(path, 'wb') as f:\n",
    "            f.write(auto.struct.pack('I', len(data.shape)))\n",
    "            f.write(auto.struct.pack(f'{len(data.shape)}I', *data.shape))\n",
    "            data.tofile(f)\n",
    "        print(f\"Wrote {path.stat().st_size :,d} bytes to {path}\")\n",
    "\n",
    "    Save(\n",
    "        index,\n",
    "        path=auto.pathlib.Path.cwd().joinpath(\n",
    "            f\"data\",\n",
    "            f\"observation\",\n",
    "            f\"OSPGeometricModel.index.vec1uc.bin\",\n",
    "        ),\n",
    "    )\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f441e25",
   "metadata": {},
   "outputs": [],
   "source": [
    "def scope():\n",
    "    # clingman's dome\n",
    "    lat0, lng0 = 35.572785, -83.518351\n",
    "    lat1, lng1 = 35.553081, -83.487097\n",
    "\n",
    "    lat = (lat0 + lat1) / 2\n",
    "    lng = (lng0 + lng1) / 2\n",
    "\n",
    "    print(f\"\"\"{lat=}, {lng=}\"\"\")\n",
    "\n",
    "    x, y, z = convert_location_to_spatial(\n",
    "        lat=lat,\n",
    "        lng=lng,\n",
    "        alt=10_000,\n",
    "    )\n",
    "    print(f'{x = :10.3f}, {y = :10.3f}, {z = :10.3f}')\n",
    "\n",
    "/scope\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce26de5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61279cf4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "gist": {
   "data": {
    "description": "share/notebooks/Sunrise - Render Preprocessed Data.ipynb",
    "public": false
   },
   "id": ""
  },
  "hide_input": false,
  "kernelspec": {
   "display_name": "venv",
   "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.10.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": false,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}