如何获取网格的纬度和经度？

Question

我目前正在开展一个项目，我需要识别美国地图上网格内的所有相交点。我的目标是将地图分成1公里的正方形部分，并获取每个正方形四个角的经纬度坐标。

假设我们有美国大陆地图，上面有网格（比如说，像经纬网）。网格每条线的边长为1公里。

如果我们放大，该网格将如下所示，点代表网格的相交点。我需要这些点的纬度和经度。

我目前的方法是使用 QGIS 生成覆盖整个美国大陆领土的 0.01 度网格。最初，我按照教程为整个世界创建多边形网格。之后，我通过从 Census.gov 导入 USA Shapefile 调整了步骤，仅关注美国。但是，我在尝试提取网格交叉点的坐标时遇到了问题。不幸的是，我一直在使用从点获取纬度、经度的方法没有返回有效的纬度和经度对。

我非常感谢任何使用此类网格点准确绘制整个美国大陆地图的建议或替代方法。谢谢！

Answer 1

找到了一种方法，通过从美国人口普查网站获取城市地区的 ESRI Shapefile。这是我最终得到的脚本，只需向其提供地理 ID 即可进行映射。您可以在我的要点上找到地理 ID 列表。

import fiona
from shapely.geometry import mapping, shape
from fiona.transform import transform_geom
from shapely import speedups
from shapely.geometry import Point
import math
import sys, os

def get_bbox_fixed(bounds):
    # Fixing bbox to start on 100 m coordinate
    minx = math.floor(bounds[0] / 100) * 100
    miny = math.floor(bounds[1] / 100) * 100
    maxx = math.ceil(bounds[2] / 100) * 100
    maxy = math.ceil(bounds[3] / 100) * 100
    # print([minx, miny, maxx, maxy])
    return [minx, miny, maxx, maxy]

def create_grid(bbox, poly, id):
    epsg3746 = "epsg:3746"  # Metric system covering US
    epsg4326 = "epsg:4326"  # Geography WGS 84 system
    points_in_grid = round(((bbox[2] + 30 - bbox[0]) / 30) * ((bbox[3] + 30 - bbox[1]) / 30))
    rows_in_grid = round((bbox[2] + 30 - bbox[0]) / 30)
    print("Number of points in BBOX: " + str(points_in_grid))
    print("Number of rows in BBOX: " + str(rows_in_grid))
    with open('output/' + id + '.csv', 'w') as out:
        row = 0
        for x in range(bbox[0], bbox[2] + 30, 30):
            for y in range(bbox[1], bbox[3] + 30, 30):
                # Creating point on row, column coordinate
                p = Point([x, y])
                # Checking if the point is inside the polygon of the urban area
                if p.within(poly):
                    # Transforming point into WGS 84
                    p_transformed = transform_geom(epsg3746, epsg4326, mapping(p))
                    p_transformed_shape = shape(p_transformed)
                    out.write(str(p_transformed_shape.x) + ',' + str(p_transformed_shape.y) + '\n')
            row += 1
            if row % 10 == 0:
                print("Processed " + str(row) + " rows from " + str(rows_in_grid) + " rows.")

speedups.enable()

epsg3746 = "epsg:3746"
epsg4269 = "epsg:4269"

if len(sys.argv) != 2:
    print("You have to specify GEOID20 of the area")
    exit()

if os.path.exists("tl_2023_us_uac20.shp"):
    with fiona.open("tl_2023_us_uac20.shp", "r") as uacs:
        # Loop via polygons of urban areas
        for uac in uacs:
            if uac['properties']['GEOID20'] == sys.argv[1]:
                print("Processing " + sys.argv[1])
                # Transforms polygon into metric system
                geom_transformed = transform_geom(epsg4269, epsg3746, uac["geometry"])
                poly = shape(geom_transformed)
                fixed_bounds = get_bbox_fixed(poly.bounds)
                create_grid(fixed_bounds, poly, uac['properties']['GEOID20'])
                break
else:
    print("You have to place us_uac_20 ESRI Shapefile in the same directory where is the script located. You can download the data from: https://www2.census.gov/geo/tiger/TIGER2023/UAC/")