pygeos-shapely-empty-geometries.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Empty geometry handling in Shapely / PyGEOS\n",
    "\n",
    "Some relevant issues:\n",
    "\n",
    "* [Constructing Empty Polygon Creates Empty Geometry Collection, or something else... · Issue #338 · Toblerity/Shapely](https://github.com/Toblerity/Shapely/issues/338 \"Constructing Empty Polygon Creates Empty Geometry Collection, or something else... · Issue #338 · Toblerity/Shapely\")\n",
    "* [Improvements to handling of empty geometries by snorfalorpagus · Pull Request #352 · Toblerity/Shapely](https://github.com/Toblerity/Shapely/pull/352 \"Improvements to handling of empty geometries by snorfalorpagus · Pull Request #352 · Toblerity/Shapely\")\n",
    "* [Empty geometries not all created the same way · Issue #742 · Toblerity/Shapely](https://github.com/Toblerity/Shapely/issues/742 \"Empty geometries not all created the same way · Issue #742 · Toblerity/Shapely\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygeos\n",
    "\n",
    "import shapely\n",
    "import shapely.geometry\n",
    "import shapely.wkt\n",
    "import shapely.wkb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shapely:  1.7.0\n",
      "PyGEOS:  0.7.1+4.gb2b845a.dirty\n",
      "GEOS:  3.8.1\n"
     ]
    }
   ],
   "source": [
    "print(\"Shapely: \", shapely.__version__)\n",
    "print(\"PyGEOS: \", pygeos.__version__)\n",
    "print(\"GEOS: \", pygeos.geos_version_string)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creation of empty geometries in Shapely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_geometry_characteristics(geom):\n",
    "    print(\"Object type      : \", type(geom).__name__)\n",
    "    try:\n",
    "        print(\"Mapping          : \",  g1.__geo_interface__)\n",
    "    except:\n",
    "        print(\"Mapping          : <ERROR>\")\n",
    "    print(\"GEOSGeometry type: \", geom.geom_type)\n",
    "    print(\"WKT              : \", geom.wkt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For polygons:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Polygon\n",
      "Mapping          :  {'type': 'Polygon', 'coordinates': ()}\n",
      "GEOSGeometry type:  GeometryCollection\n",
      "WKT              :  GEOMETRYCOLLECTION EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 1: using an empty constructor\n",
    "g1 = shapely.geometry.Polygon()\n",
    "print_geometry_characteristics(g1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Polygon\n",
      "Mapping          :  {'type': 'Polygon', 'coordinates': ()}\n",
      "GEOSGeometry type:  GeometryCollection\n",
      "WKT              :  GEOMETRYCOLLECTION EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 2: loading GeoJSON like mapping\n",
    "g2 = shapely.geometry.shape({'type': 'Polygon', 'coordinates': ()})\n",
    "print_geometry_characteristics(g2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Polygon\n",
      "Mapping          :  {'type': 'Polygon', 'coordinates': ()}\n",
      "GEOSGeometry type:  Polygon\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 3: loading empty WKB\n",
    "g3 = shapely.wkb.loads('010300000000000000', hex=True)\n",
    "print_geometry_characteristics(g3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Polygon\n",
      "Mapping          :  {'type': 'Polygon', 'coordinates': ()}\n",
      "GEOSGeometry type:  Polygon\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 4: loading empty WKT\n",
    "g4 = shapely.wkt.loads('POLYGON EMPTY')\n",
    "print_geometry_characteristics(g4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Polygon\n",
      "Mapping          :  {'type': 'Polygon', 'coordinates': ()}\n",
      "GEOSGeometry type:  Polygon\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 5: from a spatial operation (eg an empty intersection of polygons)\n",
    "g5 = shapely.geometry.box(0, 0, 1, 1).intersection(shapely.geometry.box(2, 2, 3, 3))\n",
    "print_geometry_characteristics(g5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Similarly for Points:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Point\n",
      "Mapping          : <ERROR>\n",
      "GEOSGeometry type:  GeometryCollection\n",
      "WKT              :  GEOMETRYCOLLECTION EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 1: using an empty constructor\n",
    "g1 = shapely.geometry.Point()\n",
    "print_geometry_characteristics(g1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method 2: loading GeoJSON like mapping\n",
    "# this errors on construction (IndexError)\n",
    "# g2 = shapely.geometry.shape({'type': 'Point', 'coordinates': ()})\n",
    "# print_geometry_characteristics(g2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Point\n",
      "Mapping          : <ERROR>\n",
      "GEOSGeometry type:  Point\n",
      "WKT              :  POINT EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 3: loading empty WKB\n",
    "# -> note: this is the WKB for Point(Nan, NaN), but PostGIS, GPKG, and others, use this for empty points\n",
    "#    It seems Shapely actually loads this as empty Point?\n",
    "g3 = shapely.wkb.loads('0101000000000000000000F87F000000000000F87F', hex=True)\n",
    "print_geometry_characteristics(g3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Point\n",
      "Mapping          : <ERROR>\n",
      "GEOSGeometry type:  Point\n",
      "WKT              :  POINT EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 4: loading empty WKT\n",
    "g4 = shapely.wkt.loads('POINT EMPTY')\n",
    "print_geometry_characteristics(g4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  Point\n",
      "Mapping          : <ERROR>\n",
      "GEOSGeometry type:  Point\n",
      "WKT              :  POINT EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 5: from a spatial operation (eg empty intersection of a point with a polygon)\n",
    "g5 = shapely.geometry.Point(0, 0).intersection(shapely.geometry.box(2, 2, 3, 3))\n",
    "print_geometry_characteristics(g5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Additionally, Shapely also added a \"EmptyGeometry\" class:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  EmptyGeometry\n",
      "Mapping          : <ERROR>\n",
      "GEOSGeometry type:  GeometryCollection\n",
      "WKT              :  GEOMETRYCOLLECTION EMPTY\n"
     ]
    }
   ],
   "source": [
    "g = shapely.geometry.base.EmptyGeometry()\n",
    "print_geometry_characteristics(g)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Summary**\n",
    "\n",
    "- The Python constructor always uses a GeometryCollection for empty geometries, while the other ways (from WKB/WKT, from spatial operation) creates an emtpy geometry of the actual geometry type.\n",
    "- All objects \"know\" they are a Polygon or a Point (their Python class), but this contradicts the actual GEOS geometry type when created with the Python constructor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using PyGEOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def print_geometry_characteristics(geom):\n",
    "    print(\"Object type      : \", type(geom).__name__)\n",
    "    print(\"GEOSGeometry type: \", pygeos.GeometryType(pygeos.get_type_id(geom)))\n",
    "    print(\"Is emtpy?        : \", pygeos.is_empty(geom))\n",
    "    print(\"WKT              : \", pygeos.to_wkt(geom))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For polygons:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IllegalArgumentException: point array must contain 0 or >1 elements\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Method 1: using an empty constructor\n",
    "# -> not yet supported in PyGEOS -> https://github.com/pygeos/pygeos/issues/149\n",
    "try:\n",
    "    g1 = pygeos.polygons(np.empty((0, 2)))\n",
    "    print_geometry_characteristics(g1)\n",
    "except Exception as exc:\n",
    "    print(exc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method 2: loading GeoJSON like mapping\n",
    "# Not yet implemented"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  GEOSGeometry\n",
      "GEOSGeometry type:  GeometryType.POLYGON\n",
      "Is emtpy?        :  True\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 3: loading empty WKB\n",
    "g3 = pygeos.from_wkb('010300000000000000')\n",
    "print_geometry_characteristics(g3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  GEOSGeometry\n",
      "GEOSGeometry type:  GeometryType.POLYGON\n",
      "Is emtpy?        :  True\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 4: loading empty WKT\n",
    "g4 = pygeos.from_wkt('POLYGON EMPTY')\n",
    "print_geometry_characteristics(g4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object type      :  GEOSGeometry\n",
      "GEOSGeometry type:  GeometryType.POLYGON\n",
      "Is emtpy?        :  True\n",
      "WKT              :  POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "# Method 5: from a spatial operation (eg an empty intersection of polygons)\n",
    "g5 = pygeos.intersection(pygeos.box(0, 0, 1, 1), pygeos.box(2, 2, 3, 3))\n",
    "print_geometry_characteristics(g5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The support of GEOS for \"esoteric\" empty types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The comment at https://github.com/Toblerity/Shapely/issues/742#issuecomment-508246485 shows an overview of the WKT and WKB representations of the empty version of all the geometry types in PostGIS. \n",
    "\n",
    "A quick check if GEOS supports those (at least when constructing from WKB or WKT):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pygeos.Geometry POLYGON EMPTY>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# empty polygon with third dimension\n",
    "p = pygeos.from_wkt(\"POLYGON Z EMPTY\")\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pygeos.is_empty(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pygeos.has_z(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pygeos.Geometry POLYGON EMPTY>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# empty polygon with third dimension\n",
    "p = pygeos.from_wkb(\"010300008000000000\")\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pygeos.is_empty(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pygeos.has_z(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So at least in the WKT/WKB reader, GEOS does not preserve the Z dimension for empty geometries.\n",
    "\n",
    "Let's try to create an empty polygon manually from an empty linearring with 3 dimensions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "from shapely.geos import lgeos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "cs = lgeos.GEOSCoordSeq_create(0, 3)\n",
    "lr = lgeos.GEOSGeom_createLinearRing(cs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "ob = shapely.geometry.base.geom_factory(lr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LINEARRING EMPTY\n"
     ]
    }
   ],
   "source": [
    "print(ob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ctypes import c_void_p, POINTER\n",
    "\n",
    "geos_shell = lr \n",
    "geos_holes = POINTER(c_void_p)()\n",
    "L = 0\n",
    "geos_polygon = lgeos.GEOSGeom_createPolygon(\n",
    "                c_void_p(geos_shell), geos_holes, L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = shapely.geometry.base.geom_factory(geos_polygon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "POLYGON EMPTY\n"
     ]
    }
   ],
   "source": [
    "print(poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly.has_z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppoly = pygeos.from_shapely(poly)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pygeos.has_z(ppoly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this seems to be a limitation of GEOS."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (geo-dev)",
   "language": "python",
   "name": "geo-dev"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}