GeoJSON with ipyvolume

geojson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GeoJSON on a Globe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows how to render GeoJSON data on a 3D globe. It can render coordinates and lines, but not solid polygons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geojson\n",
    "import ipyvolume as ipv\n",
    "from ipyvolume.datasets import UrlCached\n",
    "from PIL import Image\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sphere"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sphere(x, y, z, radius, num=30, color=None, texture=None, wireframe=False):\n",
    "    \"\"\"Create a sphere mesh with origin at x, y, z and radius.\n",
    "    \"\"\"\n",
    "    assert num > 0\n",
    "    u = np.linspace(0, 1, num)\n",
    "    v = np.linspace(0, 1, num)\n",
    "    u, v = np.meshgrid(u, v)\n",
    "    phi = u * 2 * np.pi\n",
    "    theta = v * np.pi\n",
    "    x = x + radius * np.cos(phi) * np.sin(theta)\n",
    "    y = y + radius * np.cos(theta)\n",
    "    z = z + radius * np.sin(phi) * np.sin(theta)\n",
    "\n",
    "    kwargs = dict(color=color or \"blue\", texture=texture, wireframe=wireframe)\n",
    "    return ipv.plot_mesh(x, y, z, u=u, v=v, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = ipv.figure()\n",
    "sphere(0, 0, 0, radius=1, num=15)\n",
    "ipv.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## GeoJSON coordinates on a globe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Render all coordinates extracted from GeoJSON data on a globe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pi_div_180 = np.pi / 180\n",
    "\n",
    "def lonlat2xyz(lon, lat, radius=1, unit='deg'):\n",
    "    \"Convert lat/lon pair to Cartesian x/y/z triple.\"\n",
    "\n",
    "    if unit == 'deg':\n",
    "        lat = lat * pi_div_180\n",
    "        lon = lon * pi_div_180\n",
    "    cos_lat = np.cos(lat)\n",
    "    x = radius * cos_lat * np.sin(lon)\n",
    "    y = radius * np.sin(lat)\n",
    "    z = radius * cos_lat * np.cos(lon)\n",
    "    return (x, y, z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "data = geojson.load(open('globe.geojson'))\n",
    "lon, lat = np.array(list(geojson.utils.coords(data))).T\n",
    "z = [lonlat2xyz(xi, yi) for (xi, yi) in zip(lon, lat)]\n",
    "x, y, z = np.array(z).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ipv.figure()\n",
    "sphere(0, 0, 0, radius=1, num=12) # , texture=image, num=100)\n",
    "ipv.scatter(x, y, z, size=1, color='limegreen', marker='sphere')\n",
    "ipv.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = \"https://www.evl.uic.edu/pape/data/Earth/1024/BigEarth.jpg\"\n",
    "image_file = UrlCached(url)\n",
    "image = Image.open(image_file.fetch())\n",
    "\n",
    "fig = ipv.figure()\n",
    "sphere(0, 0, 0, radius=1, num=24, texture=image)\n",
    "ipv.scatter(x, y, z, size=0.5, color='limegreen', marker='sphere')\n",
    "ipv.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, something is still upside down..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## GeoJSON lines on a globe\n",
    "\n",
    "Render lines extracted from GeoJSON data on a globe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def line_coords(gj):\n",
    "    \"\"\"Yield lines in given GeoJSON/TopoJSON object as lists of [lon, lat] coordinates.\n",
    "\n",
    "    Ignores Point and MultiPoint types because they don't form lines.\n",
    "    \"\"\"\n",
    "    gj_type = gj['type']\n",
    "\n",
    "    if gj_type in ['LineString']:\n",
    "        yield gj['coordinates']\n",
    "    elif gj_type in ['MultiLineString', 'Polygon']:\n",
    "        yield gj['coordinates']\n",
    "    elif gj_type == 'MultiPolygon':\n",
    "        for poly in gj['coordinates']:\n",
    "            for line in poly:\n",
    "                yield line\n",
    "    elif gj_type == 'GeometryCollection':\n",
    "        for geom in gj['geometries']:\n",
    "            for line in line_coords(geom):\n",
    "                yield line\n",
    "    elif gj_type == 'FeatureCollection':\n",
    "        for feat in gj['features']:\n",
    "            for line in line_coords(feat):\n",
    "                yield line\n",
    "    elif gj_type == 'Feature':\n",
    "        geom = gj['geometry']\n",
    "        for line in line_coords(geom):\n",
    "            yield line\n",
    "    elif gj_type == 'Topology':  # TopoJSON\n",
    "        transform = gj.get('transform', {})\n",
    "        scale = transform.get('scale', [1.0, 1.0])\n",
    "        translate = transform.get('translate', [0.0, 0.0])\n",
    "        for arc in gj['arcs']:\n",
    "            line = []\n",
    "            prev = [0, 0]\n",
    "            for point in arc:\n",
    "                prev[0] += point[0]\n",
    "                prev[1] += point[1]\n",
    "                line.append((prev[0] * scale[0] + translate[0], prev[1] * scale[1] + translate[1]))\n",
    "            yield line\n",
    "    elif gj_type in ['Point', 'MultiPoint']:\n",
    "        pass\n",
    "    else:\n",
    "        msg = f'Unknown GeoJSON/TopoJSON type: {gj_type}'\n",
    "        raise ValueError(msg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyz = []\n",
    "for line in line_coords(data):\n",
    "    if len(np.array(line).shape) < 2:\n",
    "        line = line[0]\n",
    "    L = len(line)\n",
    "    if L == 1:\n",
    "        line = line[0]\n",
    "    z = [lonlat2xyz(lon, lat) for (lon, lat) in line]\n",
    "    x, y, z = np.array(z).T\n",
    "    xyz.append([x, y, z])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "ipv.figure()\n",
    "sphere(0, 0, 0, radius=1, num=60)\n",
    "for x, y, z in xyz:\n",
    "    ipv.plot(x, y, z, color=\"limegreen\")\n",
    "ipv.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[screenshot](screenshot/geojson.gif)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}