栅格化原理
把某个点根据经纬度放在整数经纬度记录的格子里，并把格子编号与点对应起来。
第一步确定每个格子的长和宽，即经度变化量和纬度变换量：
假设测试点的经纬度是(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.864748latStart = 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 npdef 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] = Truereturn 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 distancedef 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)).kmdist_lon = distance.distance((lat_min, lon_min), (lat_min, lon_max)).kmlat_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] = Truereturn 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)
