Python地图栅格化实例 - CNPolaris
source link: https://www.cnblogs.com/cnpolaris/p/16789522.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
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()
本文作者:CNPolaris
本文链接:https://www.cnblogs.com/cnpolaris/p/16789522.html
版权声明:本作品采用知识共享署名-非商业性使用-禁止演绎 2.5 中国大陆许可协议进行许可。
Recommend
About Joyk
Aggregate valuable and interesting links.
Joyk means Joy of geeK