在使用底图创建的加利福尼亚州地图上创建等距坐标

Question

我已经使用 Basemap python 库和这个 shapefile.

创建了加利福尼亚州地图

代码如下

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(20,20))
map = Basemap(llcrnrlon=-124.48,llcrnrlat=32.53,urcrnrlon=-114.13,urcrnrlat=42.01,
             resolution='c', projection='lcc', lat_0 =  36.778259, lon_0 = -119.417)
#westlimit=-124.48; southlimit=32.53; eastlimit=-114.13; northlimit=42.01
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='#f2f2f2',lake_color='#46bcec')
map.drawcoastlines()
map.readshapefile('./CA_Counties/CA_Counties_TIGER', 'CA_Counties_TIGER')
plt.show()

输出

现在我接下来要做的是仅在加利福尼亚地图上绘制等距点，这将填满整个加利福尼亚地图（没有海洋和附近的州）。示例如下所示，因为网格中的所有点（几乎）都填满了网格。

我认为可以做到的方法是将加利福尼亚视为一个多边形或多边形，并在其中生成等间距的经度和纬度。我在看 https://gis.stackexchange.com/questions/207731/how-to-generate-random-coordinates-in-a-multipolygon-in-python 但它并没有真正解决我的问题，因为当我运行以下代码生成点数

import fiona
from shapely.geometry import shape
import random
from shapely.geometry import Point
from shapely.geometry import Polygon

def generate_random(number, polygon):
    list_of_points = []
    minx, miny, maxx, maxy = polygon.bounds
    counter = 0
    while counter < number:
        pnt = Point(random.uniform(minx, maxx), random.uniform(miny, maxy))
        if polygon.contains(pnt):
            list_of_points.append(pnt)
            counter += 1
    return list_of_points
all_points=[]
for pol in fiona.open('./CA_Counties/CA_Counties_TIGER.shp'):
    #l.append(pol['geometry'])
    all_points.append(generate_random(50, Polygon(pol['geometry']['coordinates'])))

它给了我以下错误

ValueError: A LinearRing must have at least 3 coordinate tuples

此外，我什至不确定上面的代码是否有效，是否会为我提供所有等距点（纬度和经度）以填充整个加利福尼亚地图。有没有其他方法可以做到这一点，或者有人可以用上面的代码帮助我。它必须生成纬度和经度。所以我可以在底图中用 map.plot() 函数绘制它们。还想获得数据结构中的所有点。例如列表的列表或数组的数组或任何其他运作良好的列表（可能是字典）

Answer 1

我没有完全按照您的示例进行操作，尤其是为什么您对点坐标使用随机生成。如果我理解正确的话，您希望沿着规则的网格生成点，但是只能在给定的 shapefile 中生成点。

我的建议如下：

import numpy as np
from shapely.geometry import Point
from shapely.ops import cascaded_union
import geopandas as gpd

def generate_grid_in_polygon(spacing, polygon):
    ''' This Function generates evenly spaced points within the given GeoDataFrame.
        The parameter 'spacing' defines the distance between the points in coordinate units. '''
    
    # Convert the GeoDataFrame to a single polygon
    poly_in = cascaded_union([poly for poly in polygon.geometry])

    # Get the bounds of the polygon
    minx, miny, maxx, maxy = poly_in.bounds    
    
    # Now generate the entire grid
    x_coords = list(np.arange(np.floor(minx), int(np.ceil(maxx)), spacing))
    y_coords = list(np.arange(np.floor(miny), int(np.ceil(maxy)), spacing))
    
    grid = [Point(x) for x in zip(np.meshgrid(x_coords, y_coords)[0].flatten(), np.meshgrid(x_coords, y_coords)[1].flatten())]
    
    # Finally only keep the points within the polygon
    list_of_points = [point for point in grid if point.within(poly_in)]

    return list_of_points

函数的使用示例如下：

shape_in = gpd.read_file('path/to/shapefile.shp')

points_in_poly = generate_grid_in_polygon(0.25, shape_in)
points_in_poly_gdf = gpd.GeoDataFrame(geometry=points_in_poly)

ax1 = shape_in.plot(facecolor='none', edgecolor='k')
ax1 = points_in_poly_gdf.plot(ax=ax1)

对于间距为 0.25° 的加利福尼亚州，可能如下所示：

自行调整。根据选择的间距，这可能需要很长时间才能处理。您可能希望使用空间索引来提高性能。

在使用底图创建的加利福尼亚州地图上创建等距坐标

Create Equally spaced coordinates on a California State map created with basemap

python

gis

shapefile

matplotlib-basemap

shapely