I have a function that outputs a grid of points as x and y numpy arrays for interpolation, but before I interpolate, I want to use Geopandas to perform an intersection with my research boundary (otherwise half of my interpolation points fall in the ocean).
I'm generating points like this:
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
x = np.linspace(0,100,100)
y = np.linspace(0,100,100)
x, y = np.meshgrid(x, y)
x, y = x.flatten(), y.flatten()
f, ax = plt.subplots()
plt.scatter(x, y)
plt.axis('equal')
plt.show()
Is there an efficient way to convert these numpy arrays to shapely.Point([x, y]) so they can be placed in a geopandas geodataframe?
This is my current approach:
interp_points = []
index = 0
y_list = yi.tolist()
for x in xi.tolist():
interp_points.append(Point(x,y_list[index]))
index += 1
But it seems like converting to lists and then iterating is likely not a good approach for performance, and I have approximately 160,000 points.
There is no built-in way to do this with shapely, so you need to iterate through the values yourself. For doing that, this should be a rather efficient way:
In [4]: from geopandas import GeoSeries
In [5]: s = GeoSeries(map(Point, zip(x, y)))
In [6]: s.head()
Out[6]:
0 POINT (0 0)
1 POINT (1.01010101010101 0)
2 POINT (2.02020202020202 0)
3 POINT (3.03030303030303 0)
4 POINT (4.040404040404041 0)
dtype: object
In [6]: %timeit GeoSeries(map(Point, zip(x, y)))
114 ms ± 8.45 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
(or slight alternative GeoSeries(list(zip(x, y))).map(Point))
See here for some example doing this: http://geopandas.readthedocs.io/en/latest/gallery/create_geopandas_from_pandas.html
There is some (stalled) work to include this in geopandas directly: https://github.com/geopandas/geopandas/pull/75