Two 2D arrays to coordinate sets

Spatial Digger 0
Viewed 37

I have two 2-dimensional arrays:

x
array([[0, 0, 0, 0],
       [1, 1, 1, 1],
       [2, 2, 2, 2],
       [3, 3, 3, 3],
       [4, 4, 4, 4]])
y
array([[0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3]])

What I would like to do is to generate a MultiPoint dataset in Geopandas, so we get:

points = MultiPoint([(0, 0), (0, 1), (0, 2), (0, 3), (1, 0), ... ])

How is this achieved?

geopandas python
3 Answers

This seems to do the trick, please say if there's a better way.

x = x.flatten()
y = y.flatten()

xy = list(zip(x, y))

points = MultiPoint(xy)

points_gdf = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[points])
points_gdf = points_gdf.explode(index_parts=True)
Spatial Digger0

For creating a single MultiPoint, calling shapely directly is just fine. Geopandas does have a method for going straight from vectors to a GeometryArray of points: geopandas.points_from_xy, which can be collapsed to a MultiPoint with GeoSeries.unary_union(). When working with very large lists of points, geopandas runs a bit faster, but not much.

Setting up a test with 1 million points:

import numpy as np, geopandas as gpd, shapely.geometry, pandas as pd

# generate large vectors of lat/lons
x = np.random.random(size=1_000_000) * 360 - 180
y = np.random.random(size=1_000_000) * 180 - 90
crs = 'epsg:4326'

Your example runs just a few seconds longer than the geopandas equivalent:

In [16]:  %%time
    ...:
    ...: points = shapely.geometry.MultiPoint(list(zip(x.flat, y.flat)))
    ...: points_gdf = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[points])
    ...:
    ...:
CPU times: user 15.6 s, sys: 99.1 ms, total: 15.7 s
Wall time: 15.7 s

In [17]: %%time
    ...:
    ...: points = gpd.points_from_xy(x.flat, y.flat).unary_union()
    ...: points_gdf = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[points])
    ...:
    ...:
CPU times: user 12.1 s, sys: 153 ms, total: 12.3 s
Wall time: 12.3 s

If you're not collapsing into a single geometry object, but working with large vectors of points, performing spatial operations on each of them, I'd definitely recommend the direct-to-geopandas route using points_from_xy.

Michael Delgado0

Given answer you have provided to your own question does an explode() there is no need to generate a MultiPoint. If you really want a MultiPoint you can use MultiPoint(gpd.points_from_xy(np.ravel(x), np.ravel(y)))

import numpy as np
import geopandas as gpd

x = np.array([[0, 0, 0, 0],
       [1, 1, 1, 1],
       [2, 2, 2, 2],
       [3, 3, 3, 3],
       [4, 4, 4, 4]])
y = np.array([[0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3],
       [0, 1, 2, 3]])

gpd.GeoDataFrame(geometry=gpd.points_from_xy(np.ravel(x), np.ravel(y)))
Rob Raymond0
Related