栅格化原理

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

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

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

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

地球的(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度。

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

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

#栅格化代码 import math #定义一个测试点的经纬度 testlon = 114 testlat = 22.5 #划定栅格的划分范围 lon1 = 113.75194 lon2 = 114. lat1 = 22. lat2 = 22. latStart = min(lat1, lat2); lonStart = min(lon1, lon2); #定义每个栅格大小(单位m) accuracy = 500; #计算每个栅格的经纬度增加量大小▲Lon和▲Lat deltaLon = accuracy * 360 / (2 * math.pi *  * math.cos((lat1 + lat2) * math.pi / 360)); deltaLat = accuracy * 360 / (2 * math.pi * ); #计算测试点所在栅格的经纬度编号 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., 40., 115., 117. # 栅格大小：0.05 度 grid, latitudes, longitudes = wgs84_grid(39., 40., 115., 117., 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., 40., 115., 117. # 栅格距离：1 公里 grid, latitudes, longitudes = wgs84_grid_by_distance(39., 40., 115., 117., 1) print('Grid shape:', grid.shape) print('Latitude steps:', len(latitudes)) print('Longitude steps:', len(longitudes)) print('Grid:\n', grid)