Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support GeoPandas Polygons and MultiPolygons #1285

Merged
merged 12 commits into from
Oct 16, 2023
Merged
15 changes: 13 additions & 2 deletions datashader/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,11 @@
except Exception:
dask_cudf = None

try:
import geopandas
except Exception:
geopandas = None

try:
import spatialpandas
except Exception:
Expand Down Expand Up @@ -738,8 +743,15 @@ def polygons(self, source, geometry, agg=None):
x_range = self.x_range if self.x_range is not None else (None, None)
y_range = self.y_range if self.y_range is not None else (None, None)
source = source.cx_partitions[slice(*x_range), slice(*y_range)]
glyph = PolygonGeom(geometry)
elif spatialpandas and isinstance(source, spatialpandas.GeoDataFrame):
pass
glyph = PolygonGeom(geometry)
elif geopandas and isinstance(source, geopandas.GeoDataFrame):
from .glyphs.polygon import GeopandasPolygonGeom
x_range = self.x_range if self.x_range is not None else (None, None)
y_range = self.y_range if self.y_range is not None else (None, None)
source = source.cx[slice(*x_range), slice(*y_range)]
glyph = GeopandasPolygonGeom(geometry)
else:
raise ValueError(
"source must be an instance of spatialpandas.GeoDataFrame or \n"
Expand All @@ -748,7 +760,6 @@ def polygons(self, source, geometry, agg=None):

if agg is None:
agg = any_rdn()
glyph = PolygonGeom(geometry)
return bypixel(source, self, glyph, agg)

def quadmesh(self, source, x=None, y=None, agg=None):
Expand Down
21 changes: 19 additions & 2 deletions datashader/glyphs/points.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,11 @@
cudf = None
cuda_args = None

try:
from geopandas.array import GeometryDtype as gpd_GeometryDtype
except ImportError:
gpd_GeometryDtype = type(None)

try:
import spatialpandas
except Exception:
Expand Down Expand Up @@ -72,11 +77,23 @@ def required_columns(self):
return [self.geometry]

def compute_x_bounds(self, df):
bounds = df[self.geometry].array.total_bounds_x
col = df[self.geometry]
if isinstance(col.dtype, gpd_GeometryDtype):
# geopandas
bounds = col.total_bounds[::2]
ianthomas23 marked this conversation as resolved.
Show resolved Hide resolved
else:
# spatialpandas
bounds = col.array.total_bounds_x
return self.maybe_expand_bounds(bounds)

def compute_y_bounds(self, df):
bounds = df[self.geometry].array.total_bounds_y
col = df[self.geometry]
if isinstance(col.dtype, gpd_GeometryDtype):
# geopandas
bounds = col.total_bounds[1::2]
else:
# spatialpandas
bounds = col.array.total_bounds_y
return self.maybe_expand_bounds(bounds)

@memoize
Expand Down
109 changes: 103 additions & 6 deletions datashader/glyphs/polygon.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,47 @@
except Exception:
spatialpandas = None

try:
import geopandas
import shapely
except:
geopandas = None
shapely = None


class GeopandasPolygonGeom(_GeometryLike):
@property
def geom_dtypes(self):
from geopandas.array import GeometryDtype
return (GeometryDtype,)

@memoize
def _build_extend(self, x_mapper, y_mapper, info, append, _antialias_stage_2, _antialias_stage_2_funcs):
expand_aggs_and_cols = self.expand_aggs_and_cols(append)
map_onto_pixel = _build_map_onto_pixel_for_line(x_mapper, y_mapper)
draw_polygon = _build_draw_polygon(
append, map_onto_pixel, x_mapper, y_mapper, expand_aggs_and_cols
)

perform_extend_cpu = _build_extend_geopandas_polygon_geometry(
draw_polygon, expand_aggs_and_cols
)
geom_name = self.geometry

def extend(aggs, df, vt, bounds, plot_start=True):
sx, tx, sy, ty = vt
xmin, xmax, ymin, ymax = bounds
aggs_and_cols = aggs + info(df, aggs[0].shape[:2])
geom_array = df[geom_name].array

perform_extend_cpu(
sx, tx, sy, ty,
xmin, xmax, ymin, ymax,
geom_array, *aggs_and_cols
)

return extend


class PolygonGeom(_GeometryLike):
# spatialpandas must be available if a PolygonGeom object is created.
Expand Down Expand Up @@ -52,7 +93,7 @@ def _build_draw_polygon(append, map_onto_pixel, x_mapper, y_mapper, expand_aggs_
@expand_aggs_and_cols
def draw_polygon(
i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
offsets, values, xs, ys, yincreasing, eligible,
offsets, offset_multiplier, values, xs, ys, yincreasing, eligible,
*aggs_and_cols
):
"""Draw a polygon using a winding-number scan-line algorithm
Expand All @@ -64,8 +105,8 @@ def draw_polygon(
eligible.fill(1)

# First pass, compute bounding box of polygon vertices in data coordinates
start_index = offsets[0]
stop_index = offsets[-1]
start_index = offset_multiplier*offsets[0]
stop_index = offset_multiplier*offsets[-1]
# num_edges = stop_index - start_index - 2
poly_xmin = np.min(values[start_index:stop_index:2])
poly_ymin = np.min(values[start_index + 1:stop_index:2])
Expand Down Expand Up @@ -105,8 +146,8 @@ def draw_polygon(
# Build arrays of edges in canvas coordinates
ei = 0
for j in range(len(offsets) - 1):
start = offsets[j]
stop = offsets[j + 1]
start = offset_multiplier*offsets[j]
stop = offset_multiplier*offsets[j + 1]
for k in range(start, stop - 2, 2):
x0 = values[k]
y0 = values[k + 1]
Expand Down Expand Up @@ -201,6 +242,62 @@ def draw_polygon(
return draw_polygon


def _build_extend_geopandas_polygon_geometry(
draw_polygon, expand_aggs_and_cols
):
def extend_cpu(
sx, tx, sy, ty, xmin, xmax, ymin, ymax, geometry, *aggs_and_cols
):
ragged = shapely.io.to_ragged_array(geometry)
ianthomas23 marked this conversation as resolved.
Show resolved Hide resolved
geometry_type = ragged[0]
if geometry_type not in (shapely.GeometryType.POLYGON, shapely.GeometryType.MULTIPOLYGON):
raise NotImplementedError

coords = ragged[1].ravel()
if geometry_type == shapely.GeometryType.MULTIPOLYGON:
offsets0, offsets1, offsets2 = ragged[2]
else: # POLYGON
offsets0, offsets1 = ragged[2]
offsets2 = np.arange(len(offsets1))

extend_cpu_numba(
sx, tx, sy, ty, xmin, xmax, ymin, ymax,
coords, offsets0, offsets1, offsets2, *aggs_and_cols
)

@ngjit
@expand_aggs_and_cols
def extend_cpu_numba(sx, tx, sy, ty, xmin, xmax, ymin, ymax,
coords, offsets0, offsets1, offsets2, *aggs_and_cols):
# Pre-allocate temp arrays
max_edges = 0
n_multipolygons = len(offsets2) - 1
for i in range(n_multipolygons):
polygon_inds = offsets1[offsets2[i]:offsets2[i + 1] + 1]
for j in range(len(polygon_inds) - 1):
start = offsets0[polygon_inds[j]]
stop = offsets0[polygon_inds[j + 1]]
max_edges = max(max_edges, stop - start - 1)

xs = np.full((max_edges, 2), np.nan, dtype=np.float32)
ys = np.full((max_edges, 2), np.nan, dtype=np.float32)
yincreasing = np.zeros(max_edges, dtype=np.int8)

# Initialize array indicating which edges are still eligible for processing
eligible = np.ones(max_edges, dtype=np.int8)

for i in range(n_multipolygons):
polygon_inds = offsets1[offsets2[i]:offsets2[i + 1] + 1]
for j in range(len(polygon_inds) - 1):
start = polygon_inds[j]
stop = polygon_inds[j + 1]
draw_polygon(i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
offsets0[start:stop + 1], 2,
coords, xs, ys, yincreasing, eligible, *aggs_and_cols)

return extend_cpu


def _build_extend_polygon_geometry(
draw_polygon, expand_aggs_and_cols
):
Expand Down Expand Up @@ -269,7 +366,7 @@ def extend_cpu_numba(
stop = polygon_inds[j + 1]

draw_polygon(i, sx, tx, sy, ty, xmin, xmax, ymin, ymax,
offsets2[start:stop + 1], values,
offsets2[start:stop + 1], 1, values,
xs, ys, yincreasing, eligible, *aggs_and_cols)

return extend_cpu
7 changes: 7 additions & 0 deletions datashader/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,11 @@
except Exception:
cudf = None

try:
from geopandas.array import GeometryDtype as gpd_GeometryDtype
except ImportError:
gpd_GeometryDtype = type(None)

try:
from spatialpandas.geometry import GeometryDtype
except ImportError:
Expand Down Expand Up @@ -430,6 +435,8 @@ def dshape_from_pandas_helper(col):
return datashape.Option(datashape.DateTime(tz=tz))
elif isinstance(col.dtype, (RaggedDtype, GeometryDtype)):
return col.dtype
elif gpd_GeometryDtype and isinstance(col.dtype, gpd_GeometryDtype):
return col.dtype
dshape = datashape.CType.from_numpy_dtype(col.dtype)
dshape = datashape.string if dshape == datashape.object_ else dshape
if dshape in (datashape.string, datashape.datetime_):
Expand Down
130 changes: 130 additions & 0 deletions natural_earth.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
{
"cells": [
{
"cell_type": "markdown",
"id": "715a71bf-104a-4447-96a2-e6a56147561a",
"metadata": {},
"source": [
"# Natural Earth dataset: GeoPandas and SpatialPandas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "17d8133a-60fe-4c5a-91af-784ee4e78ee2",
"metadata": {},
"outputs": [],
"source": [
"import datashader as ds\n",
"import geopandas\n",
"import geodatasets\n",
"import spatialpandas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a55e8a5e-1039-4f63-a95c-5ca2034384b6",
"metadata": {},
"outputs": [],
"source": [
"# Load data into GeoPandas and SpatialPandas GeoDataFrames\n",
"gdf = geopandas.read_file(geodatasets.get_path(\"naturalearth.land\"))\n",
"sdf = spatialpandas.GeoDataFrame(gdf) # Convert to SpatialPandas"
]
},
{
"cell_type": "markdown",
"id": "94e7ba4c-62b7-4477-aaec-929a511f8bf8",
"metadata": {},
"source": [
"## Display both outputs"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "25b7ca9d-649e-491d-9848-1d45f92b1bf1",
"metadata": {},
"outputs": [],
"source": [
"# GeoPandas\n",
"canvas = ds.Canvas(plot_height=600, plot_width=1200)\n",
"agg = canvas.polygons(source=gdf, geometry=\"geometry\", agg=ds.count())\n",
"ds.transfer_functions.shade(agg, cmap=[\"white\", \"green\"], how=\"eq_hist\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10964acc-1b8f-4411-b4ef-45be064ce32b",
"metadata": {},
"outputs": [],
"source": [
"# SpatialPandas\n",
"canvas = ds.Canvas(plot_height=600, plot_width=1200)\n",
"agg = canvas.polygons(source=sdf, geometry=\"geometry\", agg=ds.count())\n",
"ds.transfer_functions.shade(agg, cmap=[\"white\", \"green\"], how=\"eq_hist\")"
]
},
{
"cell_type": "markdown",
"id": "f4e187a1-0bd1-4b89-b64b-0b4d86567015",
"metadata": {},
"source": [
"## Time GeoPandas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7a7c5cca-9fac-4ef5-8429-cf8a8f89384f",
"metadata": {},
"outputs": [],
"source": [
"%%timeit\n",
"agg = canvas.polygons(source=gdf, geometry=\"geometry\", agg=ds.count())"
]
},
{
"cell_type": "markdown",
"id": "92c4549d-1dcd-40d7-bcd7-a60ab194fcae",
"metadata": {},
"source": [
"## Time SpatialPandas"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0e5b956e-7163-4762-9953-5353ddd661b4",
"metadata": {},
"outputs": [],
"source": [
"%%timeit\n",
"agg = canvas.polygons(source=sdf, geometry=\"geometry\", agg=ds.count())"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Loading
Loading