栅格化原理

把某个点根据经纬度放在整数经纬度记录的格子里，并把格子编号与点对应起来。

第一步确定每个格子的长和宽，即经度变化量和纬度变换量：

假设测试点的经纬度是(114度, 22.5度)

划定栅格划分的经纬度范围（大范围）为
经度范围：lon1=113.75194度，lon2=114.624187度
纬度范围：lat1=22.447837度，lat2=22.864748度
则中间点的经纬度是((lon1+lon2)/2, (lat1+lat2)/2)

那么起始点为经度和纬度取小的那个点
(latstart,lonstart)=(min(lat1,lat2),min(lon1,lon2))

规定每个栅格的大小(单位m)为500m，
下面来计算每个栅格的经度变化量和纬度变化量

假设地球是一个半径为R的圆，
那么地球的(lon1+lon2)/2度的经度圈的半径为R，绕其走一圈，路程为2* pi* R，纬度变化量为360度；

地球的(lat1+lat2)/2度的纬度圆的半径为R* cos((lat1+lat2)/2)，绕其走一圈，路程为2* pi* R* cos((lat1+lat2)/2)，经度变化量为360度。

则可以得到一组对应关系
纵向变化路程与纬变化量对应关系：2* pi* R --> 360度；
横向变化路程与经度变化量对应关系：2* pi* R* cos((lat1+lat2)/2) --> 360度。

那么纵向路程为500m时，纬度变化量为(500* 360度)/2* pi* R
横向路程为500m时，经度变化量为(500* 360度)/2* pi* R* cos((lat1+lat2)/2)

接着计算测试点(114度, 22.5度)所在栅格的经纬度编号
经度栅格编号=(114度-(起始点经度-栅格经度变化量/2))/栅格经度变化量
这里减的是(起始点经度-栅格经度变化量/2)，因为这里假设起始点是其所在栅格的右上角点，(起始点经度-栅格经度变化量/2)是起始点所在栅格中心点的经度。
纬度栅格编号同理

计算测试点所在栅格的中心点的经纬度
中心点经度 = 测试点所在栅格编号*栅格经度变化量 +起始点所在栅格中心点经度
中心点纬度同理

#栅格化代码
import math
#定义一个测试点的经纬度
testlon = 114
testlat = 22.5

#划定栅格的划分范围
lon1 = 113.75194
lon2 = 114.624187
lat1 = 22.447837
lat2 = 22.864748

latStart = min(lat1, lat2);
lonStart = min(lon1, lon2);

#定义每个栅格大小(单位m)
accuracy = 500;

#计算每个栅格的经纬度增加量大小▲Lon和▲Lat
deltaLon = accuracy * 360 / (2 * math.pi * 6371004 * math.cos((lat1 + lat2) * math.pi / 360));
deltaLat = accuracy * 360 / (2 * math.pi * 6371004);

#计算测试点所在栅格的经纬度编号
LONCOL=divmod(float(testlon) - (lonStart - deltaLon / 2) , deltaLon)[0]
LATCOL=divmod(float(testlat) - (latStart - deltaLat / 2) , deltaLat)[0]

#计算测试点所在栅格的中心点经纬度
HBLON = LONCOL*deltaLon + (lonStart - deltaLon / 2)#格子编号*格子宽+起始横坐标-半个格子宽=格子中心横坐标
HBLAT = LATCOL*deltaLat + (latStart - deltaLat / 2)

# 测试点所在栅格经纬度编号，测试点所在栅格中心点经纬度，每个栅格的经纬度变化量
LONCOL,LATCOL,HBLON,HBLAT,deltaLon,deltaLat

按照指定度数进行栅格化

import numpy as np

def wgs84_grid(lat_min, lat_max, lon_min, lon_max, grid_size):
    """将指定的 WGS84 坐标系下的区域栅格化"""
    # 计算经度和纬度的划分数量
    lat_steps = int(np.ceil((lat_max - lat_min) / grid_size))
    lon_steps = int(np.ceil((lon_max - lon_min) / grid_size))
    
    # 构造栅格网格
    latitudes = np.linspace(lat_min, lat_max, lat_steps+1)
    longitudes = np.linspace(lon_min, lon_max, lon_steps+1)
    grid = np.zeros((lat_steps, lon_steps), dtype=bool)
    
    # 遍历每个栅格，判断其是否在指定区域内
    for i in range(lat_steps):
        for j in range(lon_steps):
            lat1, lat2 = latitudes[i], latitudes[i+1]
            lon1, lon2 = longitudes[j], longitudes[j+1]
            # 判断当前栅格是否与指定区域相交
            if lat1 <= lat_max and lat2 >= lat_min and lon1 <= lon_max and lon2 >= lon_min:
                grid[i,j] = True
    return grid, latitudes, longitudes

# 测试代码
if __name__ == '__main__':
    # 北京市范围：39.442758, 40.215446, 115.420363, 117.507645
    # 栅格大小：0.05 度
    grid, latitudes, longitudes = wgs84_grid(39.442758, 40.215446, 115.420363, 117.507645, 0.05)
    print('Grid shape:', grid.shape)
    print('Latitude steps:', len(latitudes))
    print('Longitude steps:', len(longitudes))
    print('Grid:\n', grid)

按照指定距离进行栅格化

import numpy as np
from geopy import distance

def wgs84_grid_by_distance(lat_min, lat_max, lon_min, lon_max, grid_size):
    """将指定的 WGS84 坐标系下的区域栅格化"""
    # 计算经度和纬度的划分数量
    dist_lat = distance.distance((lat_min, lon_min), (lat_max, lon_min)).km
    dist_lon = distance.distance((lat_min, lon_min), (lat_min, lon_max)).km
    lat_steps = int(np.ceil(dist_lat / grid_size))
    lon_steps = int(np.ceil(dist_lon / grid_size))
    
    # 构造栅格网格
    latitudes = np.linspace(lat_min, lat_max, lat_steps+1)
    longitudes = np.linspace(lon_min, lon_max, lon_steps+1)
    grid = np.zeros((lat_steps, lon_steps), dtype=bool)
    
    # 遍历每个栅格，判断其是否在指定区域内
    for i in range(lat_steps):
        for j in range(lon_steps):
            lat1, lat2 = latitudes[i], latitudes[i+1]
            lon1, lon2 = longitudes[j], longitudes[j+1]
            # 判断当前栅格是否与指定区域相交
            if lat1 <= lat_max and lat2 >= lat_min and lon1 <= lon_max and lon2 >= lon_min:
                grid[i,j] = True
    return grid, latitudes, longitudes

# 测试代码
if __name__ == '__main__':
    # 北京市范围：39.442758, 40.215446, 115.420363, 117.507645
    # 栅格距离：1 公里
    grid, latitudes, longitudes = wgs84_grid_by_distance(39.442758, 40.215446, 115.420363, 117.507645, 1)
    print('Grid shape:', grid.shape)
    print('Latitude steps:', len(latitudes))
    print('Longitude steps:', len(longitudes))
    print('Grid:\n', grid)