Python地图栅格化实例
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()
参考链接
python-geopandas读取、创建shapefile文件、geopandas学习教程
GeoPandas官方文档
基于arcpy实现导出区域内网格中心坐标功能
gepandas安装问题
网页标题:Python地图栅格化实例
链接地址:http://azwzsj.com/article/dsogcep.html