Python地图栅格化实例
引言
shapefile是GIS中的一种非常重要的数据类型，由ESRI开发的空间数据开放格式，目前该数据格式已经成为了GIS领域的开放标准。目前绝大多数开源以及收费的GIS软件都支持该数据类型。事实上，shapefile文件指的一种文件存储的方法，实际上该种文件是由多个文件组成的。组成一个shapefile有三种文件必不可少, '.shp','.shx','.dbf'文件。
geopandas对shapefile提供了很好的读取与写出支持。geopandas库允许对几何类型进行空间操作，其中DataFrame结构相当于GIS数据中的一张属性表，使得可以直接操作矢量数据属性表，使得python中操作地理数据更加方便。本实例通过geopandas实现对地理数据的操作。
开发准备
由于geopandas库的安装需要一些前提库，因此需要先安装一些库
pip install pipwin pipwin install gdal pipwin install fiona pip install geopandas 实测以上方法可以成功在windows下安装（注：如果在Anaconda下安装geopandas更为方便）
数据准备
该数据是一段GPS扫描数据，包含经纬度。
代码实例
环境引入
import geopandas as gp import matplotlib.pyplot as plt from shapely import geometry import math GPS数据处理
lake_original_path = 'data.txt' lake_original_data = '' lake_points = [] # 读取文件 with open(lake_original_path) as f: lake_original_data = f.read() # 处理经纬度坐标 并以Point的形式添加到list中 for xy in lake_original_data.split(';'): x, _, y = xy.partition(',') x = float(x.strip()) / 100 y = float(y.strip()) / 100 lake_points.append(geometry.Point(y, x)) 创建要素
# 创建线状要素 lake_line = geometry.LineString(lake_points) # crs指定坐标系 lake_ = gp.GeoSeries(lake_line, crs='EPSG:4326') # 保存shp文件 lake_.to_file("boundary.shp", driver='ESRI Shapefile', encoding='utf-8') # 记录边界条件 用于构建栅格 x_min, y_min, x_max, y_max = lake_line.bounds[:4] # 绘图 lake_.plot() plt.show() 
构建栅格
# 栅格大小 GRID_WIDTH = 0.009 * 2 / 100 grid_rows_num = int(math.ceil((y_max - y_min) / float(GRID_WIDTH))) grid_columns_num = int(math.ceil((x_max - x_min) / float(GRID_WIDTH))) grids = [] for r in range(grid_rows_num): for c in range(grid_columns_num): grid_4coords = [] # 左上角 x_lt = x_min + c * GRID_WIDTH y_lt = y_max - r * GRID_WIDTH # 右上角 x_rt = x_lt + GRID_WIDTH y_rt = y_lt # 左下角 x_lb = x_lt y_lb = y_lt - GRID_WIDTH # 右下角 x_rb = x_rt y_rb = y_lb # 两个三角形拼接一个栅格 grid_4coords.append(geometry.Point(x_lt,y_lt)) grid_4coords.append(geometry.Point(x_rt,y_rt)) grid_4coords.append(geometry.Point(x_rb,y_rb)) grid_4coords.append(geometry.Point(x_lb,y_lb)) grid_4coords.append(geometry.Point(x_lt,y_lt)) # 创建一个网格 grids.append(geometry.LineString(grid_4coords)) grid_ = gp.GeoSeries(grids) grid_.to_file('E:\just\海韵湖智能技术实验场\data\grids.shp',driver='ESRI Shapefile', encoding='utf-8') grid_.plot() plt.show()
要素叠加
# 要素叠加 elements = [lake_line] elements += grids elements_ = gp.GeoSeries(elements) elements_.to_file('elements.shp', driver='ESRI Shapefile', encoding='utf-8') elements_.plot() plt.show()
