Extract EPSG:4326 bounding boxes from 3D Tiles dataset and plot with ipyleaflet

3D Tiles Bounds.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import numpy as np\n",
    "import pyproj\n",
    "from ipyleaflet import Map, Polygon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot Bounding Boxes of 3D Tiles Dataset\n",
    "\n",
    "Example data: https://data.kitware.com/#user/5edea5859014a6d84eabf0f0/folder/6193fc902fa25629b9d1e3eb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = '/Users/bane/Desktop/jacksonville/tileset.json'\n",
    "# f = '/Users/bane/Desktop/3d-tiles-samples-main/1.0/TilesetWithDiscreteLOD/tileset.json'\n",
    "with open(f, 'r') as f:\n",
    "    tileset_json = json.loads(f.read())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-81.67075773,  30.33012387],\n",
       "       [-81.65604684,  30.33012387],\n",
       "       [-81.65604684,  30.31681576],\n",
       "       [-81.67075773,  30.31681576],\n",
       "       [-81.67075773,  30.33012387]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def bounds_to_polygon(xmin, xmax, ymin, ymax):\n",
    "    return np.array(\n",
    "        [\n",
    "            (xmin, ymax),\n",
    "            (xmax, ymax),\n",
    "            (xmax, ymin),\n",
    "            (xmin, ymin),\n",
    "            (xmin, ymax),  # close the loop\n",
    "        ]\n",
    "    )\n",
    "    \n",
    "\n",
    "def bounding_volume_to_ll(data: dict, transform: np.ndarray = None):\n",
    "    volume = data['boundingVolume']\n",
    "    if 'transform' in data:\n",
    "        transform = np.array(data['transform']).reshape((4, 4), order='C')\n",
    "    # TODO: handle nested transforms\n",
    "    # The coordinates should be in EPSG:4978 by convention in 3D Tiles\n",
    "    # See https://github.com/CesiumGS/3d-tiles/tree/main/specification#coordinate-reference-system-crs\n",
    "    #srid = 4978\n",
    "    if 'region' in volume:\n",
    "        xmin, ymin, xmax, ymax, _, _ = volume['region']\n",
    "        # The region bounding volume specifies bounds using a geographic coordinate system (latitude, longitude, height), specifically EPSG 4979.\n",
    "        # See https://github.com/CesiumGS/3d-tiles/tree/main/specification#coordinate-reference-system-crs\n",
    "        #srid = 4979\n",
    "        # Regions are in WGS84 radians, so simply convert to degrees\n",
    "        return bounds_to_polygon(xmin, xmax, ymin, ymax) * 180. / np.pi\n",
    "    elif 'box' in volume:\n",
    "        b = volume['box']\n",
    "        # NOTE: I don't think this is right\n",
    "        center, x, y, _ = b[0:3], b[3:6], b[6:9], b[9:12]\n",
    "        xmin = min(center[0], x[0], y[0])\n",
    "        xmax = max(center[0], x[0], y[0])\n",
    "        ymin = min(center[1], x[1], y[1])\n",
    "        ymax = max(center[1], x[1], y[1])\n",
    "        coords = bounds_to_polygon(xmin, xmax, ymin, ymax)\n",
    "        coords = np.c_[coords, np.zeros(len(coords)), np.ones(len(coords))]\n",
    "         # Apply the tranformation\n",
    "        corners = (coords @ transform)[:,:-1]  # in XYZ world\n",
    "        src = pyproj.CRS(f'EPSG:4978') # 4979 -> srid\n",
    "        dest = pyproj.CRS(f'EPSG:4326')\n",
    "        transformer = pyproj.Transformer.from_proj(src, dest, always_xy=True)\n",
    "\n",
    "        x, y, _ = transformer.transform(corners.T[0], corners.T[1], corners.T[2])\n",
    "        return np.c_[x, y]\n",
    "    elif 'sphere' in volume:\n",
    "        raise NotImplementedError\n",
    "    raise ValueError(f'Bounding volume of unknown type: {volume}')\n",
    "\n",
    "\n",
    "\n",
    "bbox_rads = bounding_volume_to_ll(tileset_json['root'])\n",
    "bbox_rads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a4a396560b904360bd14559fe44d7373",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[30.330123868680694, -81.67075772657986], controls=(ZoomControl(options=['position', 'zoom_in_text'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = Map(center=bbox_rads[0][::-1].tolist(), zoom=14)\n",
    "\n",
    "polygon = Polygon(\n",
    "    locations=bbox_rads[:,::-1].tolist(),\n",
    "    color=\"green\",\n",
    ")\n",
    "\n",
    "m.add_layer(polygon)\n",
    "\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# View Children Tiles' Footprints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "transform = np.array(tileset_json['root']['transform']).reshape((4, 4), order='C')\n",
    "\n",
    "polygons = []\n",
    "for child in tileset_json['root']['children']:\n",
    "    bbox_rads = bounding_volume_to_ll(child, transform=transform)\n",
    "    polygons.append(bbox_rads[:,0:2])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9beaf09e3811470ca46a58933aa2055c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[30.329954653633806, -81.66284911437619], controls=(ZoomControl(options=['position', 'zoom_in_text'…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = Map(center=bbox_rads[0][::-1].tolist(), zoom=14)\n",
    "for p in polygons:\n",
    "    layer = Polygon(\n",
    "        locations=p[:,::-1].tolist(),\n",
    "        color=\"green\",\n",
    "    )\n",
    "\n",
    "    m.add_layer(layer)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:django] *",
   "language": "python",
   "name": "conda-env-django-py"
  },
  "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}