数据处理——研究范围栅格化

栅格化公式图

生成的栅格图

代码：

  1 import pandas as pd
  2 import numpy as np
  3 import geopandas as gpd
  4 #用geopandas包读取深圳的行政区划shapefile文件，赋值给sz变量
  5 sz = gpd.GeoDataFrame.from_file(r'data\sz.shp',encoding = 'utf-8')
  6 sz
  7 type(sz)
  8 sz['geometry'].iloc[0]
  9 #用GeoDataFrame自带的plot函数，绘制sz
 10 sz.plot()
 11 
 12 #栅格化代码
 13 import math
 14 #定义一个测试栅格划的经纬度
 15 testlon = 114
 16 testlat = 22.5
 17 
 18 #划定栅格划分范围
 19 lon1 = 113.75194
 20 lon2 = 114.624187
 21 lat1 = 22.447837
 22 lat2 = 22.864748
 23 
 24 latStart = min(lat1, lat2);
 25 lonStart = min(lon1, lon2);
 26 
 27 #定义栅格大小(单位m)
 28 accuracy = 500;
 29 
 30 #计算栅格的经纬度增加量大小▲Lon和▲Lat
 31 deltaLon = accuracy * 360 / (2 * math.pi * 6371004 * math.cos((lat1 + lat2) * math.pi / 360));
 32 deltaLat = accuracy * 360 / (2 * math.pi * 6371004);
 33 
 34 #计算栅格的经纬度编号
 35 LONCOL=divmod(float(testlon) - (lonStart - deltaLon / 2) , deltaLon)[0]
 36 LATCOL=divmod(float(testlat) - (latStart - deltaLat / 2) , deltaLat)[0]
 37 
 38 #计算栅格的中心点经纬度
 39 #HBLON = LONCOL*deltaLon + (lonStart - deltaLon / 2)#格子编号*格子宽+起始横坐标-半个格子宽=格子中心横坐标
 40 #HBLAT = LATCOL*deltaLat + (latStart - deltaLat / 2)
 41 #以下为更正，不需要减去半个格子宽
 42 HBLON = LONCOL*deltaLon + lonStart #格子编号*格子宽+起始横坐标=格子中心横坐标
 43 HBLAT = LATCOL*deltaLat + latStart 
 44 
 45 
 46 #把算好的东西print出来看看
 47 LONCOL,LATCOL,HBLON,HBLAT,deltaLon,deltaLat
 48 
 49 #该栅格的经纬度范围
 50 HBLON - deltaLon/2,HBLON + deltaLon/2,HBLAT - deltaLat/2,HBLAT + deltaLon/2
 51 
 52 from shapely.geometry import Polygon
 53 Polygon([[HBLON - deltaLon/2,HBLAT - deltaLat/2],
 54          [HBLON - deltaLon/2,HBLAT + deltaLat/2],
 55          [HBLON + deltaLon/2,HBLAT + deltaLat/2],
 56          [HBLON + deltaLon/2,HBLAT - deltaLat/2]
 57         ])
 58 
 59 import pandas as pd
 60 import numpy as np
 61 import matplotlib as mpl
 62 import matplotlib.pyplot as plt
 63 import geopandas
 64 from shapely.geometry import Point,Polygon,shape
 65 
 66 
 67 #定义空的geopandas表
 68 data = geopandas.GeoDataFrame()
 69 
 70 #定义空的list，后面循环一次就往里面加东西
 71 LONCOL = []
 72 LATCOL = []
 73 geometry = []
 74 HBLON1 = []
 75 HBLAT1 = []
 76 
 77 #计算总共要生成多少个栅格
 78 #lon方向是lonsnum个栅格
 79 lonsnum = int((lon2-lon1)/deltaLon)+1
 80 #lat方向是latsnum个栅格
 81 latsnum = int((lat2-lat1)/deltaLat)+1
 82 
 83 for i in range(lonsnum):
 84     for j in range(latsnum):
 85 
 86         HBLON = i*deltaLon + lonStart 
 87         HBLAT = j*deltaLat + latStart
 88         #把生成的数据都加入到前面定义的空list里面
 89         LONCOL.append(i)
 90         LATCOL.append(j)
 91         HBLON1.append(HBLON)
 92         HBLAT1.append(HBLAT)
 93         
 94         #生成栅格的Polygon形状
 95         #这里我们用周围的栅格推算三个顶点的位置，否则生成的栅格因为小数点取值的问题会出现小缝，无法完美覆盖
 96         HBLON_1 = (i+1)*deltaLon + lonStart
 97         HBLAT_1 = (j+1)*deltaLat + latStart 
 98         geometry.append(Polygon([
 99         (HBLON-deltaLon/2,HBLAT-deltaLat/2),
100         (HBLON_1-deltaLon/2,HBLAT-deltaLat/2),
101         (HBLON_1-deltaLon/2,HBLAT_1-deltaLat/2),
102         (HBLON-deltaLon/2,HBLAT_1-deltaLat/2)]))
103         
104 #为geopandas文件的每一列赋值为刚刚的list
105 data['LONCOL'] = LONCOL
106 data['LATCOL'] = LATCOL
107 data['HBLON'] = HBLON1
108 data['HBLAT'] = HBLAT1
109 data['geometry'] = geometry
110 
111 data.plot()
112 
113 #筛选出深圳范围的栅格
114 grid_sz = data[data.intersects(sz.unary_union)]
115 grid_sz.plot()
116 
117 # 深圳栅格的保存
118 grid_sz.to_file(r'data/grid_sz')