Sunrise - Preprocess Observation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8988a289",
   "metadata": {
    "ExecuteTime": {
     "start_time": "2023-08-16T22:51:35.412Z"
    }
   },
   "outputs": [],
   "source": [
    "# %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": "38633a1f",
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "SPECIES = \"0000143\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaf55110",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22ac97e3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.608025Z",
     "start_time": "2023-08-16T17:03:13.433756Z"
    }
   },
   "outputs": [],
   "source": [
    "%%scope\n",
    "def string_function(arr):\n",
    "    name = arr.__class__.__name__\n",
    "    if arr.size < 1024:\n",
    "        values = auto.numpy.array_repr(auto.numpy.ravel(arr))\n",
    "    else:\n",
    "        values = '...'\n",
    "        \n",
    "    shape = auto.numpy.shape(arr)\n",
    "    dtype = arr.dtype\n",
    "    return f'{name}({values}, {shape=!r}, {dtype=!r})'\n",
    "auto.numpy.set_string_function(string_function, repr=False)\n",
    "auto.numpy.set_string_function(string_function, repr=True)\n",
    "\n",
    "auto.numpy.set_printoptions(\n",
    "    threshold=4,\n",
    "    edgeitems=0,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "545d8e35",
   "metadata": {},
   "source": [
    "### Typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b501e17",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.613481Z",
     "start_time": "2023-08-16T17:03:13.610415Z"
    }
   },
   "outputs": [],
   "source": [
    "Degree = auto.typing.NewType('Degree', float)\n",
    "Radian = auto.typing.NewType('Radian', float)\n",
    "Meter = auto.typing.NewType('Meter', float)\n",
    "Kilometer = auto.typing.NewType('Kilometer', float)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "404d2bc7",
   "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": "a0dc7442",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.623809Z",
     "start_time": "2023-08-16T17:03:13.615368Z"
    }
   },
   "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_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": "9b2b2f26",
   "metadata": {},
   "source": [
    "### Tile (zoom, x, y) to location (lat, lng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd229e66",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.642501Z",
     "start_time": "2023-08-16T17:03:13.626842Z"
    }
   },
   "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",
    "    print(\"LATLNG\")\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",
    "    \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": "markdown",
   "id": "9e5aca25",
   "metadata": {},
   "source": [
    "### Dataclasses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fc4294e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.650011Z",
     "start_time": "2023-08-16T17:03:13.644468Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, order=True, frozen=True)\n",
    "class Coord:\n",
    "    \"\"\"Tile image coordinates for mapbox/etc services\"\"\"\n",
    "    z: int\n",
    "    x: int\n",
    "    y: int\n",
    "\n",
    "    @property\n",
    "    def nw(self) -> Location:\n",
    "        \"\"\"North-West lat/lng of current tile\"\"\"\n",
    "        lat, lng = tile_coordinates_to_latlng(\n",
    "            z=self.z,\n",
    "            x=self.y,\n",
    "            y=self.x,\n",
    "            # z=self.z,\n",
    "            # x=self.x,\n",
    "            # y=self.y,\n",
    "        )\n",
    "        return Location(lat=lat, lng=lng)\n",
    "\n",
    "    @property\n",
    "    def se(self) -> Location:\n",
    "        \"\"\"South-East lat/lng of current tile\"\"\"\n",
    "        lat, lng = tile_coordinates_to_latlng(\n",
    "            z=self.z,\n",
    "            x=1+self.y,\n",
    "            y=1+self.x,\n",
    "            # z=self.z,\n",
    "            # x=1+self.x,\n",
    "            # y=1+self.y,\n",
    "        )\n",
    "        return Location(lat=lat, lng=lng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fdb5aa8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.658280Z",
     "start_time": "2023-08-16T17:03:13.651767Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, frozen=True)\n",
    "class Domain:\n",
    "    \"\"\"Lower (inclusive) and upper (inclusive) bounds of latitude/longitude\"\"\"\n",
    "    lo: Coord\n",
    "    hi: Coord\n",
    "        \n",
    "    @property\n",
    "    def coords(self) -> List[Coord]:\n",
    "        \"\"\"All coordinates within lower/upper bounds\"\"\"\n",
    "        coords = []\n",
    "        for z, x, y in auto.itertools.product(\n",
    "            range(self.lo.z, 1+self.hi.z),\n",
    "            range(self.lo.x, 1+self.hi.x),\n",
    "            range(self.lo.y, 1+self.hi.y),\n",
    "        ):\n",
    "            coords.append(Coord(\n",
    "                z=z,\n",
    "                x=x,\n",
    "                y=y,\n",
    "            ))\n",
    "        \n",
    "        return coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0bdce30",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.665765Z",
     "start_time": "2023-08-16T17:03:13.660337Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, order=True, frozen=True)\n",
    "class Location:\n",
    "    \"\"\"Geographical latitude/longitude\"\"\"\n",
    "    lat: Degree\n",
    "    lng: Degree\n",
    "    alt: Optional[Meter] = auto.dataclasses.field(default=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "241550e4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.675287Z",
     "start_time": "2023-08-16T17:03:13.667775Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, order=True, frozen=True)\n",
    "class Point:\n",
    "    x: Kilometer\n",
    "    y: Kilometer\n",
    "    z: Kilometer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ddeb72b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.683984Z",
     "start_time": "2023-08-16T17:03:13.677039Z"
    }
   },
   "outputs": [],
   "source": [
    "# @dispatch\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",
    "    z: ClassVar[int] = 21\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",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1bd68ec7",
   "metadata": {},
   "source": [
    "### URL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20dd9081",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.698189Z",
     "start_time": "2023-08-16T17:03:13.688113Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass\n",
    "class TileURL:\n",
    "    scheme: str\n",
    "    netloc: str\n",
    "    prefix: str\n",
    "    z: int\n",
    "    x: int\n",
    "    y: int\n",
    "    suffix: str\n",
    "    query: Dict[str, str]\n",
    "    fragment: str\n",
    "\n",
    "    @classmethod\n",
    "    @auto.mediocreatbest.doctest\n",
    "    def decode(cls, url: str) -> Self:\n",
    "        r\"\"\"\n",
    "        \n",
    "        >>> decode(dict, 'https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg')\n",
    "        ... # doctest: +NORMALIZE_WHITESPACE\n",
    "        {'scheme': 'https',\n",
    "         'netloc': 'atlas-stg.geoplatform.gov',\n",
    "         'prefix': ['', 'styles', 'v1', 'atlas-user', 'ck58pyquo009v01p99xebegr9', 'tiles', '256'],\n",
    "         'z': 11,\n",
    "         'x': 549,\n",
    "         'y': 806,\n",
    "         'suffix': '@2x',\n",
    "         'query': [('access_token',\n",
    "                    'pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg')],\n",
    "         'fragment': ''}\n",
    "        >>> decode(dict, 'https://atlas-stg.geoplatform.gov:443/v4/atlas-user.0000214_pink/11/549/806.png?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg')\n",
    "        ... # doctest: +NORMALIZE_WHITESPACE\n",
    "        {'scheme': 'https',\n",
    "         'netloc': 'atlas-stg.geoplatform.gov:443',\n",
    "         'prefix': ['', 'v4', 'atlas-user.0000214_pink'],\n",
    "         'z': 11,\n",
    "         'x': 549,\n",
    "         'y': 806,\n",
    "         'suffix': '.png',\n",
    "         'query': [('access_token',\n",
    "                    'pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg')],\n",
    "         'fragment': ''}\n",
    "        \n",
    "        \"\"\"\n",
    "        \n",
    "        scheme, netloc, path, query, fragment = auto.urllib.parse.urlsplit(\n",
    "            url,\n",
    "        )\n",
    "        path = path.split('/')\n",
    "        *prefix, z, x, suffix = path\n",
    "        assert (match := auto.re.match(r'^(\\d+)(.*)$', suffix)) is not None, \\\n",
    "            f'{suffix = !r}'\n",
    "        y, suffix = match.groups()\n",
    "            \n",
    "        \n",
    "        z, x, y = map(int, (z, x, y))\n",
    "        query = auto.urllib.parse.parse_qsl(\n",
    "            query,\n",
    "        )\n",
    "\n",
    "        return cls(\n",
    "            scheme=scheme,\n",
    "            netloc=netloc,\n",
    "            prefix=prefix,\n",
    "            z=z,\n",
    "            x=x,\n",
    "            y=y,\n",
    "            suffix=suffix,\n",
    "            query=query,\n",
    "            fragment=fragment,\n",
    "        )\n",
    "\n",
    "    def encode(self) -> str:\n",
    "        query = auto.urllib.parse.urlencode(\n",
    "            self.query,\n",
    "            doseq=True,\n",
    "        )\n",
    "        \n",
    "        z, x, y = map(str, (self.z, self.x, self.y))\n",
    "\n",
    "        path = '/'.join([\n",
    "            *self.prefix,\n",
    "            z,\n",
    "            x,\n",
    "            y + self.suffix,\n",
    "        ])\n",
    "\n",
    "        return auto.urllib.parse.urlunsplit((\n",
    "            self.scheme,\n",
    "            self.netloc,\n",
    "            path,\n",
    "            query,\n",
    "            self.fragment,\n",
    "        ))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3786e22e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c59a74de",
   "metadata": {},
   "source": [
    "## HTTP Request Caching"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09d0137c",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85a03fb6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.711428Z",
     "start_time": "2023-08-16T17:03:13.699847Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "@run(once=True)\n",
    "def cache():\n",
    "    path = auto.pathlib.Path.cwd() / 'temp' / 'Sunrise.shelve'\n",
    "    return auto.shelve.open(str(path), 'c')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f417f9a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.718099Z",
     "start_time": "2023-08-16T17:03:13.713068Z"
    }
   },
   "outputs": [],
   "source": [
    "@run\n",
    "def fetch(cache, /):\n",
    "    def fetch(url: str, *, tqdm=None) -> bytes:\n",
    "        key = url\n",
    "        if key not in cache:\n",
    "            if tqdm is not None:\n",
    "                tqdm.set_description(f'Cache Miss: {url}')\n",
    "\n",
    "            # with auto.requests.get(url) as r:\n",
    "            #     r.raise_for_status()\n",
    "            #     cache[key] = r.content\n",
    "            try:\n",
    "                r = auto.requests.get(url)\n",
    "                r.raise_for_status()\n",
    "                content = r.content\n",
    "            except auto.requests.exceptions.HTTPError:\n",
    "                content = 404\n",
    "            \n",
    "            cache[key] = content\n",
    "            auto.time.sleep(1)\n",
    "        \n",
    "        else:\n",
    "            if tqdm is not None:\n",
    "                tqdm.set_description('Cache Hit: {url}')\n",
    "            \n",
    "            content = cache[key]\n",
    "        \n",
    "        if isinstance(content, int):\n",
    "            # raise auto.requests.exceptions.HTTPError(content)\n",
    "            return None\n",
    "        \n",
    "        return content\n",
    "\n",
    "    return fetch"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcec0562",
   "metadata": {},
   "source": [
    "## Observations\n",
    "\n",
    "We want to get the observation category, i.e. the percentage. The goal of this section is to get the tile data for a particular species for the \"pinks\" colormap and assign a category (integer, 0 to 3, inclusive) according to the nearest \"pinks\" color. The categories are: `0` for 0-25%; `1` for 25-50%; `2` for 50-75%; and `3` for 75-100%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0e2cde4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.730036Z",
     "start_time": "2023-08-16T17:03:13.719760Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, frozen=True)\n",
    "class Observation:\n",
    "    common_name: str\n",
    "    irma_id: str\n",
    "    latin_name_ref: str\n",
    "    latin_name: str\n",
    "\n",
    "@run\n",
    "def observations(fetch, /):\n",
    "    observations = []\n",
    "    for observation in (\n",
    "        Observation(**kwargs)\n",
    "        for kwargs in (\n",
    "            auto.json.loads(fetch('https://species.atlasofthesmokies.org/atbirecords/lexical_index.json'))\n",
    "        )['items']\n",
    "    ):\n",
    "        observations.append(observation)\n",
    "    \n",
    "    return observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26d3e8c6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.737911Z",
     "start_time": "2023-08-16T17:03:13.732181Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(after=print)\n",
    "def observation(observations, /):\n",
    "    for observation in observations:\n",
    "        if observation.irma_id not in [\n",
    "            '0000223',\n",
    "        ]:\n",
    "            continue\n",
    "        \n",
    "        return observation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2989ed92",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.895641Z",
     "start_time": "2023-08-16T17:03:13.740045Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='observation')\n",
    "def db(observation, /):\n",
    "    @auto.mediocreatbest.contextgenerator(call=lambda generator: lambda sql, /, **parameters: generator.send((sql, parameters)))\n",
    "    def db(observation: Observation=observation, /):\n",
    "        root = auto.pathlib.Path('/mnt/seenas2/data/johns_nps_run/maxent/mbtiles/out')\n",
    "        path = root / f'atlas-user.{observation.irma_id}_pink.mbtiles'\n",
    "        \n",
    "        with auto.sqlite3.connect(path) as db:\n",
    "            db.row_factory = auto.sqlite3.Row\n",
    "            \n",
    "            cursor = None\n",
    "            while True:\n",
    "                sql, parameters = yield cursor\n",
    "                cursor = db.execute(sql, parameters)\n",
    "    \n",
    "    return db\n",
    "\n",
    "@run(scope='observation', after=print)\n",
    "def z():\n",
    "    # z = 15\n",
    "    z = 13\n",
    "    return z\n",
    "\n",
    "@auto.dataclasses.dataclass(eq=True, frozen=True)\n",
    "class ObservationTile:\n",
    "    coord: Coord\n",
    "    image: bytes = auto.dataclasses.field(repr=False)\n",
    "\n",
    "@run(scope='observation', after=lambda tiles: print(f'{len(tiles)=!r}'))\n",
    "def tiles(db, z, /):\n",
    "    with db() as db:\n",
    "        tiles = []\n",
    "        for tile in db(r'''\n",
    "            SELECT\n",
    "                tile_row AS y_tms,\n",
    "                tile_column AS x,\n",
    "                zoom_level AS z,\n",
    "                tile_data AS image\n",
    "            FROM tiles\n",
    "            WHERE zoom_level = :z\n",
    "        ''', z=z):\n",
    "            tiles.append(ObservationTile(\n",
    "                coord=Coord(\n",
    "                    z=tile['z'],\n",
    "                    y=2 ** tile['z'] - 1 - tile['y_tms'],\n",
    "                    x=tile['x'],\n",
    "                ),\n",
    "                image=tile['image'],\n",
    "            ))\n",
    "    \n",
    "    return tiles\n",
    "\n",
    "@run(scope='observation', after=print)\n",
    "def xs(tiles, /):\n",
    "    return sorted(set(tile.coord.x for tile in tiles))\n",
    "\n",
    "@run(scope='observation', after=print)\n",
    "def ys(tiles, /):\n",
    "    return sorted(set(tile.coord.y for tile in tiles))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b2c3662",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:13.903594Z",
     "start_time": "2023-08-16T17:03:13.897856Z"
    }
   },
   "outputs": [],
   "source": [
    "def __stitch_tiles(tiles, /, size: int=128, flip=False):\n",
    "    auto.PIL; from PIL import Image\n",
    "    \n",
    "    xs = set(tile.coord.x for tile in tiles)\n",
    "    ys = set(tile.coord.y for tile in tiles)\n",
    "    \n",
    "    xs = sorted(xs)\n",
    "    ys = sorted(ys)\n",
    "    \n",
    "    nx = len(xs)\n",
    "    ny = len(ys)\n",
    "    \n",
    "    sz = size\n",
    "    \n",
    "    composite = Image.new('RGB', (nx*sz, ny*sz))\n",
    "    \n",
    "    for tile in tiles:\n",
    "        xi = xs.index(tile.coord.x)\n",
    "        yi = ys.index(tile.coord.y)\n",
    "        \n",
    "        with auto.io.BytesIO(tile.image) as f:\n",
    "            try:\n",
    "                image = Image.open(f)\n",
    "            except auto.PIL.UnidentifiedImageError:\n",
    "                continue\n",
    "\n",
    "            image.thumbnail((sz, sz))\n",
    "\n",
    "        # composite.paste(image, (yi*sz, xi*sz))\n",
    "        if flip==True:\n",
    "            composite.paste(image, (yi*sz, xi*sz))\n",
    "        else:\n",
    "            composite.paste(image, (xi*sz, yi*sz))\n",
    "            \n",
    "    \n",
    "    return composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a275135f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:19.009430Z",
     "start_time": "2023-08-16T17:03:13.905578Z"
    }
   },
   "outputs": [],
   "source": [
    "__stitch_tiles(auto.mediocreatbest.g['observation__tiles'], size=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fed9a65a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:19.018444Z",
     "start_time": "2023-08-16T17:03:19.012292Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='observation', cond=False)\n",
    "def _(tiles, /):\n",
    "    \"\"\"Remove redundant pixels from tile data\"\"\"\n",
    "    \n",
    "    def Image(tile: Tile, /) -> auto.PIL.Image.Image:\n",
    "        image: bytes = tile.image\n",
    "        with auto.io.BytesIO(image) as f:\n",
    "            image = auto.PIL.Image.open(f)\n",
    "            image.load()\n",
    "        \n",
    "        return image\n",
    "    \n",
    "    tile = tiles[9]\n",
    "    image = Image(tile)\n",
    "    \n",
    "    image = auto.numpy.array(image)\n",
    "    assert len(image.shape) == 3 and image.shape[-1] == 4, \\\n",
    "        f\"\"\"Wrong shape: {image.shape=!r}\"\"\"\n",
    "    NR, NC, _ = image.shape\n",
    "\n",
    "    image = image[:NR, :NC, 0]\n",
    "    \n",
    "    auto.IPython.display.display(image)\n",
    "    \n",
    "    raise NotImplementedError"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e289a888",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.344698Z",
     "start_time": "2023-08-16T17:03:19.020671Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='observation')\n",
    "def interpolator(tiles, z, xs, ys, /):\n",
    "    res = 256\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",
    "    @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",
    "    @auto.mediocreatbest.coroutine  # (call=lambda gen: lambda *, z, x, y: gen.send((z, x, y)))\n",
    "    def category() -> Generator[Optional[auto.numpy.ndarray], tuple[float, float, float], None]:\n",
    "        category = None\n",
    "        while True:\n",
    "            coord = yield category\n",
    "            \n",
    "            for tile in tiles:\n",
    "                if tile.coord == coord:\n",
    "                    break\n",
    "            \n",
    "            else:\n",
    "#                 print(f'No tile found: {coord=!r}')\n",
    "                category = None\n",
    "                continue\n",
    "\n",
    "#             (exists,) ,= db(r'''\n",
    "#                 SELECT EXISTS(SELECT 1 FROM tiles WHERE zoom_level = :z AND tile_column = :x AND tile_row = :y)\n",
    "#             ''', z=z, x=x, y=y)\n",
    "\n",
    "#             if not exists:\n",
    "#                 category = None\n",
    "#                 continue\n",
    "\n",
    "#             (image,) ,= db(r'''\n",
    "#                 SELECT tile_data FROM tiles WHERE zoom_level = :z AND tile_column = :x AND tile_row = :y\n",
    "#             ''', z=z, x=x, y=y)\n",
    "\n",
    "            with auto.io.BytesIO(tile.image) as f:\n",
    "                image = auto.PIL.Image.open(f)\n",
    "                image.load()\n",
    "            \n",
    "            image = image.convert('RGB')\n",
    "            image = yield image\n",
    "\n",
    "            rgb = auto.numpy.array(image)\n",
    "\n",
    "            assert rgb.shape == (res, res, 3), \\\n",
    "                f\"\"\"Wrong shape for observation tile image: {rgb.shape=!r}\"\"\"\n",
    "\n",
    "            rgb = rgb.reshape((res*res, 3))\n",
    "\n",
    "            diffs = auto.numpy.subtract(rgb[:, None, :], pinks[None, :, :], dtype=auto.numpy.float32)\n",
    "            assert diffs.shape == (res*res, len(pinks), 3)\n",
    "\n",
    "            dists = auto.numpy.linalg.norm(diffs, axis=-1)\n",
    "            assert dists.shape == (res*res, len(pinks))\n",
    "\n",
    "            category = auto.numpy.argmin(dists, axis=-1)\n",
    "            assert category.shape == (res*res,), \\\n",
    "                f\"\"\"Wrong shape: {category.shape=!r}\"\"\"\n",
    "\n",
    "            category = category.astype(auto.numpy.uint8)\n",
    "            category = category.reshape((res, res))\n",
    "\n",
    "    @auto.mediocreatbest.coroutine\n",
    "    def debounce(coroutine, /, *, rate: 'item/seconds'):\n",
    "        def now() -> 'seconds':\n",
    "            return auto.time.clock_gettime(auto.time.CLOCK_MONOTONIC)\n",
    "        \n",
    "        prev = float('-inf')\n",
    "        while True:\n",
    "            value = yield\n",
    "            if (curr := now()) - prev > 1.0 / rate:\n",
    "                prev = curr\n",
    "                coroutine.send(value)\n",
    "    \n",
    "    dz = int(auto.math.log2(res))\n",
    "    assert z + dz == Tixel.z\n",
    "    points = (\n",
    "        auto.numpy.arange(\n",
    "            (min(ys) + 0) * 2 ** dz,\n",
    "            (max(ys) + 1) * 2 ** dz,\n",
    "        ),\n",
    "        auto.numpy.arange(\n",
    "            (min(xs) + 0) * 2 ** dz,\n",
    "            (max(xs) + 1) * 2 ** dz,\n",
    "        ),\n",
    "    )\n",
    "    print(f'{points[0][[0, -1]]=!r}')\n",
    "    print(f'{points[1][[0, -1]]=!r}')\n",
    "    \n",
    "    values = auto.numpy.zeros(\n",
    "        shape=(\n",
    "            len(ys)*res,\n",
    "            len(xs)*res,\n",
    "        ),\n",
    "        dtype=auto.numpy.uint8,\n",
    "    )\n",
    "    \n",
    "    with auto.contextlib.ExitStack() as stack:\n",
    "        category = stack.enter_context( auto.contextlib.closing(category()) )\n",
    "        display = stack.enter_context( auto.contextlib.closing(display()) )\n",
    "        display = stack.enter_context(auto.contextlib.closing( debounce(display, rate=1) ))\n",
    "        \n",
    "        for (yi, y), (xi, x) in auto.tqdm.tqdm(auto.itertools.product(\n",
    "            enumerate(ys),\n",
    "            enumerate(xs),\n",
    "        ), total=len(ys)*len(xs)):\n",
    "            image = category.send(Coord(\n",
    "                z=z,\n",
    "                x=x,\n",
    "                y=y,\n",
    "            ))\n",
    "\n",
    "            if image is None:\n",
    "                continue\n",
    "            \n",
    "            display.send(image)\n",
    "            \n",
    "            src = category.send(image)\n",
    "                \n",
    "            values[\n",
    "                slice(\n",
    "                    (yi + 0) * res,\n",
    "                    (yi + 1) * res,\n",
    "                ),\n",
    "                slice(\n",
    "                    (xi + 0) * res,\n",
    "                    (xi + 1) * res,\n",
    "                ),\n",
    "            ] = src\n",
    "    \n",
    "    return auto.scipy.interpolate.RegularGridInterpolator(\n",
    "        points=points,\n",
    "        values=values,\n",
    "        method='nearest',\n",
    "        bounds_error=False,\n",
    "        fill_value=0,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b7666df",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.353542Z",
     "start_time": "2023-08-16T17:03:51.347943Z"
    }
   },
   "outputs": [],
   "source": [
    "auto.mediocreatbest.g['observation__interpolator'].grid[0][[0, -1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a56d5578",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.375469Z",
     "start_time": "2023-08-16T17:03:51.355668Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='observation')\n",
    "def __view(interpolator, /):\n",
    "    rs, cs = interpolator.grid\n",
    "    \n",
    "    def decimate(\n",
    "        arr: auto.numpy.ndarray,\n",
    "        /,\n",
    "        *,\n",
    "        length: int,\n",
    "    ) -> auto.numpy.ndarray:\n",
    "        print(f'{len(arr)=!r}')\n",
    "        stride = int(auto.math.floor( len(arr) / length ))\n",
    "        arr = arr[::stride]\n",
    "        arr = arr[:length]\n",
    "        return arr\n",
    "    \n",
    "    rs = decimate(rs, length=256)\n",
    "    cs = decimate(cs, length=256)\n",
    "    \n",
    "    grid = auto.numpy.hstack(tuple(\n",
    "        v.reshape((-1, 1))\n",
    "        for v in auto.numpy.meshgrid(\n",
    "            rs,\n",
    "            cs,\n",
    "            indexing='ij',\n",
    "        )\n",
    "    ))\n",
    "    assert grid.shape == (len(rs) * len(cs), 2)\n",
    "    \n",
    "    confidence = (interpolator(grid, method='linear') * 255.0 / 4.0).astype(auto.numpy.uint8)\n",
    "    confidence = confidence.reshape((len(rs), len(cs)))\n",
    "    \n",
    "    image = auto.PIL.Image.fromarray(confidence, mode='L')\n",
    "    \n",
    "    auto.IPython.display.display(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f0d0c93",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.396720Z",
     "start_time": "2023-08-16T17:03:51.377473Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "## @Tile.from_.register(Location)\n",
    "@auto.mediocreatbest.doctest\n",
    "def convert_location_to_tile(*, lat: Degree, lng: Degree, zoom: int, math=auto.math) -> tuple[int, int]:\n",
    "    r\"\"\"\n",
    "    \n",
    "    >>> convert_location_to_tile(lat=35.6, lng=-83.52, zoom=11)\n",
    "    (807, 548)\n",
    "    >>> convert_location_to_tile(lat=35.0, lng=-84.5, zoom=15)\n",
    "    (12979, 8692)\n",
    "    >>> convert_location_to_tile(lat=36.25, lng=-82.5, zoom=15)\n",
    "    (12839, 8874)\n",
    "    \n",
    "    \"\"\"\n",
    "    # Thanks https://stackoverflow.com/a/72476578\n",
    "    \n",
    "    #> lat_rad = math.radians(lat_deg)\n",
    "    lat_rad = math.radians(lat)\n",
    "    \n",
    "    #> n = 2.0 ** zoom\n",
    "    n = 2.0 ** zoom\n",
    "    \n",
    "    #> xtile = int((lon_deg + 180.0) / 360.0 * n)\n",
    "    col = int((lng + 180.0) / 360.0 * n)\n",
    "    \n",
    "    #> ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)\n",
    "    row = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)\n",
    "    \n",
    "    #> return (xtile, ytile)\n",
    "    return row, col\n",
    "\n",
    "\n",
    "@run(scope='observation')\n",
    "def __view(interpolator, /):\n",
    "    lat, lng = (\n",
    "        (35.6201452, -83.4372325)  # Clingman's Dome\n",
    "    )\n",
    "    \n",
    "    r, c = convert_location_to_tile(\n",
    "        lat=lat,\n",
    "        lng=lng,\n",
    "        zoom=Tixel.z,\n",
    "    )\n",
    "    print(f'{r=!r} {c=!r}')\n",
    "    \n",
    "    k = 2\n",
    "    rs = auto.numpy.arange(r-k*128, r+k*128+1, k)\n",
    "    cs = auto.numpy.arange(c-k*128, c+k*128+1, k)\n",
    "    \n",
    "#     rs, cs = interpolator.grid\n",
    "    \n",
    "#     def decimate(\n",
    "#         arr: auto.numpy.ndarray,\n",
    "#         /,\n",
    "#         *,\n",
    "#         length: int,\n",
    "#     ) -> auto.numpy.ndarray:\n",
    "#         print(f'{len(arr)=!r}')\n",
    "#         stride = int(auto.math.floor( len(arr) / length ))\n",
    "#         arr = arr[::stride]\n",
    "#         arr = arr[:length]\n",
    "#         return arr\n",
    "    \n",
    "#     rs = decimate(rs, length=256)\n",
    "#     cs = decimate(cs, length=256)\n",
    "    \n",
    "    grid = auto.numpy.hstack(tuple(\n",
    "        v.reshape((-1, 1))\n",
    "        for v in auto.numpy.meshgrid(\n",
    "            rs,\n",
    "            cs,\n",
    "            indexing='ij',\n",
    "        )\n",
    "    ))\n",
    "    assert grid.shape == (len(rs) * len(cs), 2)\n",
    "    \n",
    "    confidence = (interpolator(grid, method='linear').astype(auto.numpy.float32) * 255.0 / 4.0).astype(auto.numpy.uint8)\n",
    "    confidence = confidence.reshape((len(rs), len(cs)))\n",
    "    \n",
    "    image = auto.PIL.Image.fromarray(confidence, mode='L')\n",
    "    \n",
    "    auto.IPython.display.display(image)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7a1caa2",
   "metadata": {},
   "source": [
    "## Background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d499464",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.512973Z",
     "start_time": "2023-08-16T17:03:51.398971Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(after=auto.pprint.pp, scope='background')\n",
    "def domain():\n",
    "    from bs4 import BeautifulSoup\n",
    "    \n",
    "    images = ''\n",
    "    template = '<img alt=\"\" role=\"presentation\" src=\"https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{x}/{y}\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 512px; height: 512px; transform: translate3d(439px, -19px, 0px); opacity: 1;\">'\n",
    "    zoom = 5\n",
    "    # for x in range((2 ** zoom)-1):\n",
    "    #     for y in range((2 ** zoom)-10):\n",
    "    #         imgurl = template.format(z = zoom, x = x, y = y)\n",
    "    #         images = images + imgurl\n",
    "    #         print(imgurl)\n",
    "    # soup = BeautifulSoup(images, 'html.parser')\n",
    "    \n",
    "    soup = BeautifulSoup(r'''\n",
    "        <img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, 237px, 0px); opacity: 1;\">\n",
    "        <img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, -275px, 0px); opacity: 1;\">\n",
    "        <img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, -19px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/807@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, 237px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/805@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, -275px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/548/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(439px, 493px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(695px, 493px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/547/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(183px, 493px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/550/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(951px, 493px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/546/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(-73px, 493px, 0px); opacity: 1;\"><img alt=\"\" role=\"presentation\" src=\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/551/808@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\" class=\"leaflet-tile leaflet-tile-loaded\" style=\"width: 256px; height: 256px; transform: translate3d(1207px, 493px, 0px); opacity: 1;\">\n",
    "    ''', 'html.parser')\n",
    "    \n",
    "    \n",
    "    urls = []\n",
    "    for img in soup('img'):\n",
    "        src = img['src']\n",
    "        url = TileURL.decode(src)\n",
    "        urls.append(url)\n",
    "\n",
    "    return Domain(\n",
    "        lo=Coord(z=11, x=546-1, y=805-1),\n",
    "        hi=Coord(z=11, x=551+1, y=808+1),\n",
    "    )\n",
    "    \n",
    "    # return Domain(\n",
    "    #     lo=Coord(\n",
    "    #         z=min(url.z for url in urls),\n",
    "    #         x=min(url.x for url in urls),\n",
    "    #         y=min(url.y for url in urls),\n",
    "    #     ),\n",
    "    #     hi=Coord(\n",
    "    #         z=max(url.z for url in urls),\n",
    "    #         x=max(url.x for url in urls),\n",
    "    #         y=max(url.y for url in urls),\n",
    "    #     ),\n",
    "    # )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f497d51",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.547423Z",
     "start_time": "2023-08-16T17:03:51.514805Z"
    }
   },
   "outputs": [],
   "source": [
    "@auto.dataclasses.dataclass(eq=True, frozen=True)\n",
    "class Tile:\n",
    "    coord: Coord\n",
    "    url: str = auto.dataclasses.field(repr=False)\n",
    "    image: bytes = auto.dataclasses.field(repr=False)\n",
    "\n",
    "@run(scope='background', after=lambda x: print(len(x)))\n",
    "def tiles(domain, fetch, /):\n",
    "    def make_url(coord: Coord, /, *, url=TileURL.decode(\n",
    "        # r\"\"\"https://atlas-stg.geoplatform.gov/styles/v1/atlas-user/ck58pyquo009v01p99xebegr9/tiles/256/11/549/806@2x?access_token=pk.eyJ1IjoiYXRsYXMtdXNlciIsImEiOiJjazFmdGx2bjQwMDAwMG5wZmYwbmJwbmE2In0.lWXK2UexpXuyVitesLdwUg\"\"\"\n",
    "        r\"\"\"https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/0/0/0\"\"\"\n",
    "    )):\n",
    "        return auto.dataclasses.replace(\n",
    "            url,\n",
    "            z=coord.z,\n",
    "            # x=coord.x,\n",
    "            # y=coord.y,\n",
    "            x=coord.y,\n",
    "            y=coord.x,\n",
    "        )\n",
    "    \n",
    "    tiles = []\n",
    "    for coord in (tqdm := auto.tqdm.tqdm(domain.coords)):\n",
    "        url = make_url(coord).encode()\n",
    "        # print(f'url: {url}')\n",
    "        # print(f'coord: {coord.x}, {coord.y}')\n",
    "        # print(f'url: {url}')\n",
    "        image = fetch(url, tqdm=tqdm)\n",
    "\n",
    "        tiles.append(Tile(\n",
    "            coord=coord,\n",
    "            url=url,\n",
    "            image=image,\n",
    "        ))\n",
    "\n",
    "    return tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5495a0e6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.557769Z",
     "start_time": "2023-08-16T17:03:51.555441Z"
    }
   },
   "outputs": [],
   "source": [
    "__stitch_tiles(run.g['background__tiles'], size=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9765ff6d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.701897Z",
     "start_time": "2023-08-16T17:03:51.560168Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='background')\n",
    "def interpolator(tiles, /):\n",
    "    # res = 512\n",
    "    res = 256\n",
    "    \n",
    "    zs = sorted(set(tile.coord.z for tile in tiles))\n",
    "    assert len(zs) == 1\n",
    "    ys = sorted(set(tile.coord.y for tile in tiles))\n",
    "    xs = sorted(set(tile.coord.x for tile in tiles))\n",
    "    \n",
    "    print(f\"\"\"{len(ys)=!r} {len(xs)=!r}\"\"\")\n",
    "    \n",
    "    # This is cursed.\n",
    "    points = (\n",
    "        2 ** (Tixel.z - auto.math.log2(res) - zs[0]) * (auto.numpy.repeat(res * auto.numpy.array(ys)[:, None], res, axis=1) + auto.numpy.arange(res)[None, :]).ravel(),\n",
    "        2 ** (Tixel.z - auto.math.log2(res) - zs[0]) * (auto.numpy.repeat(res * auto.numpy.array(xs)[:, None], res, axis=1) + auto.numpy.arange(res)[None, :]).ravel(),\n",
    "    )\n",
    "    \n",
    "    print(f\"\"\"ys: {points[0][[0, -1]]=!r}\n",
    "xs: {points[1][[0, -1]]=!r}\"\"\")\n",
    "    \n",
    "    values = auto.numpy.zeros(\n",
    "        shape=(len(points[0]), len(points[1]), 3),\n",
    "        dtype=auto.numpy.float32,\n",
    "    )\n",
    "    print(f\"\"\"{values.shape=!r} {values.dtype=!r}\"\"\")\n",
    "\n",
    "    for tile in tiles:\n",
    "        yi = ys.index(tile.coord.y)\n",
    "        xi = xs.index(tile.coord.x)\n",
    "        \n",
    "        with auto.io.BytesIO(tile.image) as f:\n",
    "            image = auto.PIL.Image.open(f)\n",
    "            image.load()\n",
    "            image = auto.PIL.ImageOps.flip(image)\n",
    "            \n",
    "        image = image.convert('RGB')\n",
    "        image = auto.numpy.array(image, dtype=auto.numpy.float32)\n",
    "        image = image[::-1,:,:]\n",
    "        assert image.shape == (res, res, 3), image.shape\n",
    "        assert image.dtype == auto.numpy.float32\n",
    "        assert 1.0 < image.max() <= 255.0\n",
    "        \n",
    "        image /= 255.0\n",
    "        assert 0.0 < image.max() <= 1.0\n",
    "        \n",
    "#         import numpy.lib.recfunctions\n",
    "#         image = auto.numpy.lib.recfunctions.unstructured_to_structured(\n",
    "#             image,\n",
    "#             dtype=values.dtype,\n",
    "#         )\n",
    "#         assert image.shape == (res, res), \\\n",
    "#             f\"\"\"Wrong shape: {image.shape=!r}\"\"\"\n",
    "#         assert image.dtype == values.dtype\n",
    "        \n",
    "        values[\n",
    "            slice(\n",
    "                (yi+0)*res,\n",
    "                (yi+1)*res,\n",
    "            ),\n",
    "            slice(\n",
    "                (xi+0)*res,\n",
    "                (xi+1)*res,\n",
    "            ),\n",
    "        ] = image\n",
    "               \n",
    "    \n",
    "    return auto.scipy.interpolate.RegularGridInterpolator(\n",
    "        points,\n",
    "        values,\n",
    "        bounds_error=False,\n",
    "        fill_value=0,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8ea267f",
   "metadata": {},
   "source": [
    "## Elevation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afb50050",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:51.779075Z",
     "start_time": "2023-08-16T17:03:51.704783Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='elevation', after=print)\n",
    "# @auto.mediocreatbest.doctest\n",
    "def domain(z=5):\n",
    "    # z = 2\n",
    "    # r'''\n",
    "    # >>> domain(z=0)\n",
    "    # Domain(lo=Coord(z=0, x=0, y=0), hi=Coord(z=0, x=0, y=0))\n",
    "    # >>> domain(z=1)\n",
    "    # Domain(lo=Coord(z=1, x=0, y=0), hi=Coord(z=1, x=1, y=1))\n",
    "    # >>> domain(z=2)\n",
    "    # Domain(lo=Coord(z=2, x=0, y=0), hi=Coord(z=2, x=3, y=3))\n",
    "    # >>> sum(1 for _ in domain(z=0).coords)\n",
    "    # 1\n",
    "    # >>> sum(1 for _ in domain(z=1).coords)\n",
    "    # 4\n",
    "    # >>> sum(1 for _ in domain(z=2).coords)\n",
    "    # 16\n",
    "    # >>> sum(1 for _ in domain(z=3).coords)\n",
    "    # 64\n",
    "    # '''\n",
    "    return Domain(\n",
    "        lo=Coord(z=11, x=546-1, y=805-1),\n",
    "        hi=Coord(z=11, x=551+1, y=808+1),\n",
    "        # lo=Coord(z=z, x=0, y=0),\n",
    "        # hi=Coord(z=z, x=(2**z)-1, y=(2**z)-10),\n",
    "    )\n",
    "\n",
    "@run(scope='elevation')\n",
    "def tiles(domain, fetch, /):\n",
    "    def make_url(coord: Coord, /, *, url=TileURL.decode(\n",
    "        r\"\"\"https://api.mapbox.com/v4/mapbox.terrain-rgb/14/3826/6127@2x.pngraw?access_token=pk.eyJ1IjoidGhvYnNvbjIiLCJhIjoiY2oxZmdqbnQzMDBpbjJxczR6dWoyemUxNiJ9.SEBHSdHLP_lZGD43r-_IDQ\"\"\"\n",
    "    )) -> URL:\n",
    "        return auto.dataclasses.replace(\n",
    "            url,\n",
    "            z=coord.z,\n",
    "            x=coord.x,\n",
    "            y=coord.y,\n",
    "        )\n",
    "    \n",
    "    tiles = []\n",
    "    for coord in (tqdm := auto.tqdm.tqdm(domain.coords)):\n",
    "        url = make_url(coord).encode()\n",
    "        image = fetch(url, tqdm=tqdm)\n",
    "        \n",
    "        tiles.append(Tile(\n",
    "            coord=coord,\n",
    "            url=url,\n",
    "            image=image,\n",
    "        ))\n",
    "    \n",
    "    return tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d110898b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:52.358733Z",
     "start_time": "2023-08-16T17:03:51.796226Z"
    }
   },
   "outputs": [],
   "source": [
    "__stitch_tiles(auto.mediocreatbest.g['elevation__tiles'], size=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6b2b479",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:52.905833Z",
     "start_time": "2023-08-16T17:03:52.361164Z"
    }
   },
   "outputs": [],
   "source": [
    "@run(scope='elevation')\n",
    "def interpolator(tiles, /):\n",
    "    res = 512\n",
    "    \n",
    "    def tile_altitude(tile: Tile, /) -> auto.numpy.NDArray[auto.numpy.float32]:\n",
    "        with auto.io.BytesIO(tile.image) as f:\n",
    "            image = auto.PIL.Image.open(f)\n",
    "            image = auto.PIL.ImageOps.flip(image) # flip the tile vertically because we reverse the order as well\n",
    "            image.load()\n",
    "        \n",
    "        image = image.convert('RGB')\n",
    "        \n",
    "        rgb = auto.numpy.array(image)\n",
    "        assert len(rgb.shape) == 3, \\\n",
    "            f\"\"\"Expected dimension 3: {rgb.shape=!r}\"\"\"\n",
    "        assert rgb.shape[0] == rgb.shape[1], \\\n",
    "            f\"\"\"Expected square image: {rgb.shape=!r}\"\"\"\n",
    "        assert rgb.shape[0] == res, \\\n",
    "            f\"\"\"Expected resolution: {res}; Got resolution: {rgb.shape[0]}\"\"\"\n",
    "        assert rgb.shape[2] == 3, \\\n",
    "            f\"\"\"Expected last component to be rgb values: {rgb.shape=!r}\"\"\"\n",
    "        assert rgb.dtype == auto.numpy.uint8, \\\n",
    "            f\"\"\"Expected dtype to be uint8_t: {rgb.dtype=!r}\"\"\"\n",
    "        \n",
    "        # Thanks https://docs.mapbox.com/data/tilesets/guides/access-elevation-data/#decode-data\n",
    "        #> elevation = -10000 + (({R} * 256 * 256 + {G} * 256 + {B}) * 0.1)\n",
    "        \n",
    "        altitude = auto.numpy.zeros(\n",
    "            shape=(rgb.shape[0], rgb.shape[1]),\n",
    "            dtype=auto.numpy.float32,\n",
    "        )\n",
    "        \n",
    "        altitude += rgb[:, :, 0] * (1 * 256 * 256)\n",
    "        altitude += rgb[:, :, 1] * (1 * 256)\n",
    "        altitude += rgb[:, :, 2] * (1)\n",
    "        altitude *= 0.1\n",
    "        altitude -= 10000\n",
    "        \n",
    "        return altitude\n",
    "\n",
    "    xs = sorted(set(tile.coord.x for tile in tiles))\n",
    "    # Reverse sort elevation tiles on the y because for some reason they orient upside down from typical mercator proj \n",
    "    ys = sorted(set(tile.coord.y for tile in tiles))\n",
    "    ys = ys[::-1]\n",
    "    zs = sorted(set(tile.coord.z for tile in tiles))\n",
    "    \n",
    "    assert len(zs) == 1\n",
    "    z ,= zs\n",
    "    \n",
    "    print(f'{z=!r}')\n",
    "    \n",
    "    dz = int(auto.math.log2(res))\n",
    "    print(f'{dz=!r} {res=!r}')\n",
    "    points = (\n",
    "        auto.numpy.arange(\n",
    "            (min(ys) + 0) * 2 ** dz,\n",
    "            (max(ys) + 1) * 2 ** dz,\n",
    "        ),\n",
    "        auto.numpy.arange(\n",
    "            (min(xs) + 0) * 2 ** dz,\n",
    "            (max(xs) + 1) * 2 ** dz,\n",
    "        ),\n",
    "    )\n",
    "    z += dz\n",
    "    print(f'{z=!r}')\n",
    "    \n",
    "    print(f'{points[0][0]=!r}')\n",
    "    print(f'{points[0][-1]=!r}')\n",
    "    \n",
    "    if z < Tixel.z:\n",
    "        dz = Tixel.z - z\n",
    "        print(f'{dz=!r}')\n",
    "        points = (\n",
    "            points[0] * 2 ** dz,\n",
    "            points[1] * 2 ** dz,\n",
    "        )\n",
    "        z += dz\n",
    "        print(f'{z=!r}')\n",
    "        \n",
    "        print(f'{points[0][0]=!r}')\n",
    "        print(f'{points[0][-1]=!r}')\n",
    "    \n",
    "    assert z == Tixel.z\n",
    "    \n",
    "    values = auto.numpy.zeros(\n",
    "        shape=(\n",
    "            len(ys)*res,\n",
    "            len(xs)*res,\n",
    "        ),\n",
    "        dtype=auto.numpy.float32,\n",
    "    )\n",
    "    \n",
    "    for tile in tiles:\n",
    "        xi = xs.index(tile.coord.x)\n",
    "        yi = ys.index(tile.coord.y)\n",
    "        \n",
    "        values[\n",
    "            slice(\n",
    "                (yi + 0) * res,\n",
    "                (yi + 1) * res,\n",
    "            ),\n",
    "            slice(\n",
    "                (xi + 0) * res,\n",
    "                (xi + 1) * res,\n",
    "            ),\n",
    "        ] = tile_altitude(tile)\n",
    "\n",
    "    return auto.scipy.interpolate.RegularGridInterpolator(\n",
    "        points=points,\n",
    "        values=values,\n",
    "        bounds_error=False,\n",
    "        fill_value=0,\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bfc6b15",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:03:53.846206Z",
     "start_time": "2023-08-16T17:03:52.908242Z"
    }
   },
   "outputs": [],
   "source": [
    "auto.mediocreatbest.g['elevation__interpolator']([\n",
    "    auto.mediocreatbest.g['elevation__interpolator'].grid[0][100],\n",
    "    auto.mediocreatbest.g['elevation__interpolator'].grid[1][100],\n",
    "], method='quintic')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e24c44f",
   "metadata": {},
   "source": [
    "## Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01449241",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:28:08.202738Z",
     "start_time": "2023-08-16T17:28:06.107846Z"
    }
   },
   "outputs": [],
   "source": [
    "@run\n",
    "def __visualize(elevation__interpolator, observation__interpolator, background__interpolator, /):\n",
    "    def visualize(\n",
    "        *,\n",
    "        name: str,\n",
    "        interpolator: auto.scipy.interpolate.RegularGridInterpolator,\n",
    "        num: int,\n",
    "    ):\n",
    "        lat0, lng0 = 35.572785, -83.518351\n",
    "        lat1, lng1 = 35.553081, -83.487097\n",
    "        \n",
    "        r0, c0 = convert_location_to_tile(\n",
    "            lat=lat0,\n",
    "            lng=lng0,\n",
    "            zoom=Tixel.z,\n",
    "        )\n",
    "        \n",
    "        r1, c1 = convert_location_to_tile(\n",
    "            lat=lat1,\n",
    "            lng=lng1,\n",
    "            zoom=Tixel.z,\n",
    "        )\n",
    "        \n",
    "        rlo, rhi = min(r0, r1), max(r0, r1)\n",
    "        clo, chi = min(c0, c1), max(c0, c1)\n",
    "\n",
    "        # real_rlo = interpolator.grid[0].min()\n",
    "        # real_rhi = interpolator.grid[0].max()\n",
    "        # real_clo = interpolator.grid[1].min()\n",
    "        # real_chi = interpolator.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",
    "\n",
    "        \n",
    "        rs = auto.numpy.linspace(rlo, rhi, num=num)\n",
    "        cs = auto.numpy.linspace(clo, chi, num=num)\n",
    "        \n",
    "        grid = auto.numpy.hstack(tuple(\n",
    "            v.reshape((-1, 1))\n",
    "            for v in auto.numpy.meshgrid(\n",
    "                rs,\n",
    "                cs,\n",
    "                indexing='ij',\n",
    "            )\n",
    "        ))\n",
    "        assert grid.shape == (len(cs)*len(rs), 2), \\\n",
    "            f\"\"\"Wrong shape: {grid.shape=!r}\"\"\"\n",
    "        \n",
    "        values = interpolator(grid)\n",
    "        values = values.reshape((len(rs), len(cs), -1))\n",
    "        \n",
    "        %matplotlib inline\n",
    "        plt = auto.matplotlib.pyplot\n",
    "        \n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(111)\n",
    "        im = ax.imshow(values)\n",
    "        ax.set_yticks(auto.numpy.arange(0, len(rs), 5))\n",
    "        ax.set_yticklabels([f'{int(r)}' for r in rs[::5]], rotation=0)\n",
    "        ax.set_xticks(auto.numpy.arange(0, len(cs), 5))\n",
    "        ax.set_xticklabels([f'{int(c)}' for c in cs[::5]], rotation=90)\n",
    "        \n",
    "        ax.set_title(f\"Sample of {name} around Clingman's Dome\")\n",
    "        ax.set_xlabel('Tixel Column Number')\n",
    "        ax.set_ylabel('Tixel Row Number', rotation=90)\n",
    "        \n",
    "        fig.colorbar(im)\n",
    "        fig.tight_layout()\n",
    "        plt.show()\n",
    "    \n",
    "    visualize(\n",
    "        name='elevation',\n",
    "        interpolator=elevation__interpolator,\n",
    "        num=128,\n",
    "    )\n",
    "    \n",
    "    visualize(\n",
    "        name='observation',\n",
    "        interpolator=observation__interpolator,\n",
    "        num=128,\n",
    "    )\n",
    "    \n",
    "    visualize(\n",
    "        name='background',\n",
    "        interpolator=background__interpolator,\n",
    "        num=128,\n",
    "    )\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae88466e",
   "metadata": {},
   "source": [
    "### Save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b22e54fd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-16T17:04:15.273285Z",
     "start_time": "2023-08-16T17:03:53.848170Z"
    }
   },
   "outputs": [],
   "source": [
    "@run\n",
    "def __save(elevation__interpolator, observation__interpolator, background__interpolator, /):\n",
    "    def save(\n",
    "        *,\n",
    "        name: str,\n",
    "        interpolator: auto.scipy.interpolate.RegularGridInterpolator,\n",
    "    ):\n",
    "        path = auto.pathlib.Path.cwd() / 'temp' / f'park-{name}.npz'\n",
    "        auto.numpy.savez_compressed(\n",
    "            path,\n",
    "            points0=interpolator.grid[0],\n",
    "            points1=interpolator.grid[1],\n",
    "            values=interpolator.values,\n",
    "        )\n",
    "        \n",
    "        print(f'Wrote {path.stat().st_size :,d} bytes to {path}')\n",
    "    \n",
    "    save(\n",
    "        name='elevation',\n",
    "        interpolator=elevation__interpolator,\n",
    "    )\n",
    "    \n",
    "    save(\n",
    "        name='observation',\n",
    "        interpolator=observation__interpolator,\n",
    "    )\n",
    "\n",
    "    \n",
    "    save(\n",
    "        name='background',\n",
    "        interpolator=background__interpolator,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80ed2294",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "gist": {
   "data": {
    "description": "share/notebooks/Sunrise.ipynb",
    "public": false
   },
   "id": ""
  },
  "hide_input": false,
  "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.10.12"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "539px",
    "width": "751px"
   },
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "265.4px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}