Skip to content

Instantly share code, notes, and snippets.

@jorisvandenbossche

jorisvandenbossche/spherely-geopandas-integration.ipynb

Created February 26, 2025 22:34
Show Gist options
  • Save jorisvandenbossche/2e8ecee0d109b77e03c45fca7bac385a to your computer and use it in GitHub Desktop.
Save jorisvandenbossche/2e8ecee0d109b77e03c45fca7bac385a to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
spherely-geopandas-integration.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "d6c5c1b3-b91b-4f89-bce9-317f69ee3f70",
"metadata": {},
"source": [
"# Showcasing in-progress GeoPandas - spherely integration"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "445adf13-4144-4e87-8f0d-0ccf129d6a03",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Running cmake --build & --install in /home/joris/scipy/repos/spherely/build/skbuild2\n"
]
}
],
"source": [
"import geodatasets\n",
"import geopandas\n",
"import shapely\n",
"import spherely"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "e5dd9125-87ac-4edd-979f-08a4f13c19f8",
"metadata": {},
"outputs": [],
"source": [
"land = geopandas.read_file(geodatasets.get_path(\"naturalearth.land\"))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b51d6958-97ef-4717-bf11-d6c87af24bce",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>featurecla</th>\n",
" <th>scalerank</th>\n",
" <th>min_zoom</th>\n",
" <th>geometry</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-59.57209 -80.04018, -59.86585 -80.5...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-159.20818 -79.49706, -161.1276 -79....</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>0.0</td>\n",
" <td>POLYGON ((-45.15476 -78.04707, -43.92083 -78.4...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-121.21151 -73.50099, -119.91885 -73...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-125.55957 -73.48135, -124.03188 -73...</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" featurecla scalerank min_zoom \\\n",
"0 Land 1 1.0 \n",
"1 Land 1 1.0 \n",
"2 Land 1 0.0 \n",
"3 Land 1 1.0 \n",
"4 Land 1 1.0 \n",
"\n",
" geometry \n",
"0 POLYGON ((-59.57209 -80.04018, -59.86585 -80.5... \n",
"1 POLYGON ((-159.20818 -79.49706, -161.1276 -79.... \n",
"2 POLYGON ((-45.15476 -78.04707, -43.92083 -78.4... \n",
"3 POLYGON ((-121.21151 -73.50099, -119.91885 -73... \n",
"4 POLYGON ((-125.55957 -73.48135, -124.03188 -73... "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land.head()"
]
},
{
"cell_type": "markdown",
"id": "95e8b4bc-0366-42c8-88f9-8a66ee104f59",
"metadata": {},
"source": [
"Converting to a _spherical_ geometry dtype:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "93a5688a-fe9f-4a9f-9abc-1bd5d84b9cce",
"metadata": {},
"outputs": [],
"source": [
"land_geog = land.astype({\"geometry\": geopandas.array.GeometryDtype(engine=\"spherical\")})"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "05cb5b9d-6efd-40ee-95ff-74db14b5e281",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>featurecla</th>\n",
" <th>scalerank</th>\n",
" <th>min_zoom</th>\n",
" <th>geometry</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-59.5721 -80.04018, -60.09111 -79.83...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-160.16789 -79.56563, -159.20818 -79...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>0.0</td>\n",
" <td>POLYGON ((-44.53779 -78.26259, -45.15476 -78.0...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-120.56518 -73.57936, -121.21151 -73...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>Land</td>\n",
" <td>1</td>\n",
" <td>1.0</td>\n",
" <td>POLYGON ((-124.79572 -73.67731, -125.55957 -73...</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" featurecla scalerank min_zoom \\\n",
"0 Land 1 1.0 \n",
"1 Land 1 1.0 \n",
"2 Land 1 0.0 \n",
"3 Land 1 1.0 \n",
"4 Land 1 1.0 \n",
"\n",
" geometry \n",
"0 POLYGON ((-59.5721 -80.04018, -60.09111 -79.83... \n",
"1 POLYGON ((-160.16789 -79.56563, -159.20818 -79... \n",
"2 POLYGON ((-44.53779 -78.26259, -45.15476 -78.0... \n",
"3 POLYGON ((-120.56518 -73.57936, -121.21151 -73... \n",
"4 POLYGON ((-124.79572 -73.67731, -125.55957 -73... "
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.head()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "d017cdb7-a1c2-439d-9c95-99b6d2050e2c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"geometry[spherical]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.geometry.dtype"
]
},
{
"cell_type": "markdown",
"id": "192762b3-8666-4aa2-a6b2-feabdef6075d",
"metadata": {},
"source": [
"Now calculating the area works without warning and gives you the result in m²!"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ffd256f6-52a0-4b77-b212-75b1ed5bbfb6",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/tmp/ipykernel_818321/3342291871.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'area' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.\n",
"\n",
" land.geometry.area\n"
]
},
{
"data": {
"text/plain": [
"0 4.202195\n",
"1 3.716703\n",
"2 20.382879\n",
"3 1.541898\n",
"4 0.683017\n",
" ... \n",
"122 2.437153\n",
"123 13.689594\n",
"124 19.395611\n",
"125 104.819482\n",
"126 677.509565\n",
"Length: 127, dtype: float64"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land.geometry.area"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "550710fc-13da-4bda-9800-72045ad931a2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 8.588523e+09\n",
"1 8.813227e+09\n",
"2 4.552045e+10\n",
"3 5.340779e+09\n",
"4 2.386070e+09\n",
" ... \n",
"122 4.970397e+09\n",
"123 2.952446e+10\n",
"124 4.287863e+10\n",
"125 2.162656e+11\n",
"126 2.190031e+12\n",
"Length: 127, dtype: float64"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.geometry.area"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "31ef1358-1839-4acc-8cff-2adbcae4e839",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 9.443957e+06\n",
"1 1.109759e+07\n",
"2 9.073724e+06\n",
"3 1.087813e+07\n",
"4 1.100197e+07\n",
" ... \n",
"122 9.260411e+06\n",
"123 1.002959e+07\n",
"124 9.921514e+06\n",
"125 9.608350e+06\n",
"126 7.645315e+06\n",
"Length: 127, dtype: float64"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.distance(spherely.create_point(0, 0))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "8857375c-e9a7-42e9-9dc5-b09a54377929",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 POINT (-62.28401 -80.48564)\n",
"1 POINT (-161.46469 -78.93375)\n",
"2 POINT (-47.7858 -79.54124)\n",
"3 POINT (-120.92793 -73.73881)\n",
"4 POINT (-125.897 -73.56808)\n",
" ... \n",
"122 POINT (48.43643 80.51494)\n",
"123 POINT (96.26788 79.92518)\n",
"124 POINT (-91.43212 79.65842)\n",
"125 POINT (-79.75826 80.07678)\n",
"126 POINT (-41.96454 73.14658)\n",
"Length: 127, dtype: geometry"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.centroid"
]
},
{
"cell_type": "markdown",
"id": "f881c68c-51d4-4089-8c5e-9505669e3e18",
"metadata": {},
"source": [
"Methods that are not yet implemented raise an error (just need to add an informative message saying they need to cast to GeometryDtype to be able to use this method):"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "b6a148e0-a612-4beb-959b-3662a4fda9df",
"metadata": {},
"outputs": [
{
"ename": "NotImplementedError",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[11], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mland_geog\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mbuffer\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m10\u001b[39;49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m~/scipy/repos/geopandas/geopandas/base.py:5062\u001b[0m, in \u001b[0;36mGeoPandasBase.buffer\u001b[0;34m(self, distance, resolution, cap_style, join_style, mitre_limit, single_sided, **kwargs)\u001b[0m\n\u001b[1;32m 4988\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mbuffer\u001b[39m(\n\u001b[1;32m 4989\u001b[0m \u001b[38;5;28mself\u001b[39m,\n\u001b[1;32m 4990\u001b[0m distance,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 4996\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs,\n\u001b[1;32m 4997\u001b[0m ):\n\u001b[1;32m 4998\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Returns a ``GeoSeries`` of geometries representing all points within\u001b[39;00m\n\u001b[1;32m 4999\u001b[0m \u001b[38;5;124;03m a given ``distance`` of each geometric object.\u001b[39;00m\n\u001b[1;32m 5000\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 5060\u001b[0m \n\u001b[1;32m 5061\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 5062\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_delegate_geo_method\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 5063\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mbuffer\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5064\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5065\u001b[0m \u001b[43m \u001b[49m\u001b[43mdistance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdistance\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5066\u001b[0m \u001b[43m \u001b[49m\u001b[43mresolution\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mresolution\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5067\u001b[0m \u001b[43m \u001b[49m\u001b[43mcap_style\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcap_style\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5068\u001b[0m \u001b[43m \u001b[49m\u001b[43mjoin_style\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mjoin_style\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5069\u001b[0m \u001b[43m \u001b[49m\u001b[43mmitre_limit\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmitre_limit\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5070\u001b[0m \u001b[43m \u001b[49m\u001b[43msingle_sided\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msingle_sided\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5071\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 5072\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n",
"File \u001b[0;32m~/scipy/repos/geopandas/geopandas/base.py:119\u001b[0m, in \u001b[0;36m_delegate_geo_method\u001b[0;34m(op, this, **kwargs)\u001b[0m\n\u001b[1;32m 116\u001b[0m kwargs[key] \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39masarray(val)\n\u001b[1;32m 118\u001b[0m a_this \u001b[38;5;241m=\u001b[39m this\u001b[38;5;241m.\u001b[39mgeometry\u001b[38;5;241m.\u001b[39mvalues\n\u001b[0;32m--> 119\u001b[0m data \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mgetattr\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43ma_this\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mop\u001b[49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 120\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m GeoSeries(data, index\u001b[38;5;241m=\u001b[39mthis\u001b[38;5;241m.\u001b[39mindex, crs\u001b[38;5;241m=\u001b[39mthis\u001b[38;5;241m.\u001b[39mcrs)\n",
"File \u001b[0;32m~/scipy/repos/geopandas/geopandas/_array_spherical.py:703\u001b[0m, in \u001b[0;36mGeographyArray.buffer\u001b[0;34m(self, distance, resolution, **kwargs)\u001b[0m\n\u001b[1;32m 702\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mbuffer\u001b[39m(\u001b[38;5;28mself\u001b[39m, distance, resolution\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m16\u001b[39m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m--> 703\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mNotImplementedError\u001b[39;00m\n",
"\u001b[0;31mNotImplementedError\u001b[0m: "
]
}
],
"source": [
"land_geog.buffer(10)"
]
},
{
"cell_type": "markdown",
"id": "1da741fc-a939-4211-a678-69eb7245af4b",
"metadata": {},
"source": [
"Binary predicates and overlays also already work:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "227ffb51-e7ce-4d27-bf90-d102a325e9c4",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 False\n",
"1 False\n",
"2 False\n",
"3 False\n",
"4 False\n",
" ... \n",
"122 False\n",
"123 False\n",
"124 False\n",
"125 False\n",
"126 False\n",
"Length: 127, dtype: bool"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.intersects(land_geog.geometry[112])"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "2dca0506-f5da-4d4b-8e5b-f4160121195c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0 GEOMETRYCOLLECTION EMPTY\n",
"1 GEOMETRYCOLLECTION EMPTY\n",
"2 GEOMETRYCOLLECTION EMPTY\n",
"3 GEOMETRYCOLLECTION EMPTY\n",
"4 GEOMETRYCOLLECTION EMPTY\n",
" ... \n",
"122 GEOMETRYCOLLECTION EMPTY\n",
"123 GEOMETRYCOLLECTION EMPTY\n",
"124 GEOMETRYCOLLECTION EMPTY\n",
"125 GEOMETRYCOLLECTION EMPTY\n",
"126 GEOMETRYCOLLECTION EMPTY\n",
"Length: 127, dtype: geometry"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.intersection(land_geog.geometry[112])"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "ee2e53ae-55a7-4631-900a-b88963fc4bc5",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Polygon 127\n",
"Name: count, dtype: int64"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"land_geog.geometry.geom_type.value_counts()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5abc135c-bc86-42ea-88e9-9678a2facb3b",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python (spherely-dev)",
"language": "python",
"name": "spherely-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.12.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment