Skip to main content

用于处理矢量数据的地理处理实用程序

项目描述

Vector IO - 用于处理矢量数据的地理处理实用程序。

  • 蟒蛇> = 3.6
  • gdal >= 2.2
  • rar
  • 解压

描述

这个项目是一个基于GDAL处理矢量数据的工具。这个工具是一个关于 gdal 的信封,旨在快速、智能、简单地处理不同类型的矢量数据。vectorIO 提供对(读取和写入)geojson、wkt、Shapefile 和 KML 的支持,支持不同空间数据类型之间的快速切换,并为来自 gdal 的警告提供异常处理程序。

安装

码头工人

创建图像并实例化容器:

# access the directory where is the Dockerfile
docker image build -t vectorio-env:001 . # build the image
# vectorio-env:001 - can be any name with the version of the your preference
docker container run -it vectorio-env:001 # instantiate a new container

Ubuntu 18.04

  • rar
apt-get install rar unrar
  • 格达尔

在 ubuntu 上安装 gdal

  • 蟒蛇的Gdal
gdalinfo --version
pip3 install gdal==<gdal_version>

特征

读写 Geojson

使用 geojson 数据。默认情况下,数据源创建为 WGS84。

  • 准备数据
from vectorio.vector import Geojson
data = '{"type": "FeatureCollection","features": [{"type": "Feature","properties": {},"geometry": {"type": "Polygon","coordinates": [[[-44.89013671875,-6.577303118123875],[-46.29638671874999,-7.460517719883772],[-44.4287109375,-7.318881730366743],[-44.89013671875,-6.577303118123875]]]}}]}'
gjs = Geojson(data)
  • 读取所有数据
# Features
gjs.feature_collection()

# Geometries
gjs.geometry_collection()
  • 阅读和迭代每个功能
for feat in gjs.features():
    print(feat)
  • 创建一个新的 geojson 文件
gjs.write('data.geojson')
  • 从geojson文件中读取
with open('data.geojson') as f:
    gj= Geojson(f.read())
    gj.feature_collection()

读写 WKT

使用 wkt 数据。支持几何集合和单一几何。默认情况下,数据源创建为 WGS84。

wkt 对象有一些参数:

WKT(as_geometry_collection=True, srid=4326)
  • as_geometry_collection:当数据是单个几何时,返回一个几何集合,方法是集合

  • srid:WKT 的初始 SRID。

  • 准备数据

from vectorio.vector import WKT
data = "GEOMETRYCOLLECTION(POINT(-48.740641051554974 -9.249606262178954), LINESTRING(-50.278726989054974 -11.023166202413554,-48.608805114054974 -10.375450023701761))"
wkt = WKT(data)
  • 读取所有数据
wkt.geometry_collection()
  • 读取和迭代每个几何图形
for geom in wkt.geometries():
    print(geom)
  • 创建一个新的 wkt 文件
wkt.write('data.wkt')
  • 从 wkt 文件中读取
with open('data.wkt') as f:
    wkt = WKT(f.read())
    wkt.geometry_collection()

读写Shapefile

使用读写 shapefile。支持压缩为 .zip 和 .rar 的 shapefile。默认情况下,数据源是基于 .prj 文件上的投影创建的。obs:.rar 文件的读写仅适用于 linux 操作系统。只有vectorio.compress.Rar(压缩引擎)类有这个限制。其他类可用于任何操作系统。

  • 准备数据
from vectorio.vector import Shapefile
shape = Shapefile('data.shp')
  • 从 .shp 文件中读取所有数据
shape.feature_collection()
  • 从 .shp 文件中读取和迭代每个功能。
# Interanting over features
for feat in shape.features():
    print(feat)

# Interanting over geometries
for geom in shape.geometries():
    print(geom)
  • 创建一个新的 shapefile(创建文件 .shp、.shx、.dbf、.prj)
shape.write('out.shp')
# >>> out.shp
读写压缩的Shapefile

默认情况下,该算法将反复搜索压缩文件中的 .shp、.shx、.dbf、.prj 文件。该算法将搜索每个扩展名的第一个文件,如果压缩文件包含2个(或更多).shp文件,或2个(或更多).prj文件,将获得第一个.shp文件和第一个.prj文件.

  • 从 zip 处理
from vectorio.vector import Shapefile, ShapefileCompressed
from vectorio.compress import Zip

shape = ShapefileCompressed(Shapefile('data.zip'), compress_engine=Zip())
shape.feature_collection() # reading all data

for feat in shape.features():  # iterating over each item
    print(feat)

shape.write('out.zip') # Creating a shapefile compressed as .zip
# >>> out.zip
  • 从 .rar 处理(仅适用于 linux 操作系统
from vectorio.vector import Shapefile, ShapefileCompressed
from vectorio.compress import Rar

shape = ShapefileCompressed(Shapefile('data.rar'), compress_engine=Rar())
shape.feature_collection()  # read all data

for feat in shape.features():  # iterating over each item
    print(feat)

shape.write('out.rar') # Creating a shapefile compressed as .rar
# >>> out.rar

读写 KML

  • 目前,KML 作为 Geojson 处理。
 vectorio.vector 导入 KML 
kml  =  KML 'data.kml' 

# Interanting over features
for feat in shape.features():
    print(feat)

# Interanting over geometries
for geom in shape.geometries():
    print(geom)

公里'out.geojson' 

重新投影向量

空间重投影适用于实现接口 IVectorIO 的相同地理类型。如果输入 srid (in_srid) 被省略,将使用几何中的 srid。

  • 重新投影 shapefile
from vectorio.vector import Shapefile, ShapefileCompressed, VectorReprojected
from vectorio.compress import Zip
shape = VectorReprojected(
    ShapefileCompressed(Shapefile('data_utm22.zip'), compress_engine=Zip()),
    in_srid=31982, out_srid=4674
)

shape.feature_collection()  # reading all data

for feat in shape.features():  # iterating by each feature
    print(feat)

shape.write('data_reprojected.zip')  # creating a new shapefile
  • 重新投影 WKT

默认情况下,wkt 在 WGS84 空间参考中。

from vectorio.vector import WKT, VectorReprojected
data = 'POLYGON((-49.698036566343376 -9.951372897703846,-51.148231878843376 -11.591810720955946,-48.467567816343376 -11.763953408065282,-49.698036566343376 -9.951372897703846))'
wkt = VectorReprojected(WKT(data), out_srid=31982)

wkt.geometry_collection()  # reading all data

for geom in wkt.geometries():  # iterating by each geometry
    print(geom)

wkt.write('data-reprojected.wkt')  # creating a new wkt file
  • 重新投影 Geojson

默认情况下,geojson 位于 WGS84 空间参考中。

from vectorio.vector import Geojson, VectorReprojected
data = '{"type": "FeatureCollection","features": [{"type": "Feature","properties": {},"geometry": {"type": "Polygon","coordinates": [[[-45.992889404296875,-9.654907854199012],[-46.12884521484374,-9.72259300616733],[-45.96954345703125,-9.738835407948073],[-45.992889404296875,-9.654907854199012]]]}}]}'
gjs = VectorReprojected(Geojson(data), out_srid=31982)

gjs.feature_collection()  # reading all data

for feat in gjs.features():  # iterating by each feature
    print(feat)

gjs.write('data-reprojected.geojson')  # creating a new geojson file

地理数据之间的快速切换

要执行快速切换,必须使用包vectorio.vector上的VectorComposite

VectorComposite(input_vector_obj, ouput_vector_obj)
从 geojson 快速切换到 wkt
  • 准备数据
from vectorio.vector import Geojson, WKT, VectorComposite
data = '{"type": "FeatureCollection","features": [{"type": "Feature","properties": {},"geometry": {"type": "Polygon","coordinates": [[[-44.89013671875,-6.577303118123875],[-46.29638671874999,-7.460517719883772],[-44.4287109375,-7.318881730366743],[-44.89013671875,-6.577303118123875]]]}}]}'
vector = VectorComposite(Geojson(data), WKT())
  • 从 geojson 读取所有几何图形作为 wkt
vector.geometry_collection()
  • 作为 wkt 迭代所有几何图形
for geom_wkt in vector.geometries():
    print(geom_wkt)
  • 创建 wkt 文件
vector.write('output.wkt')
从 wkt 快速切换到 shapefile 作为 zip
from vectorio.vector import Shapefile, ShapefileCompressed, WKT, VectorComposite
from vectorio.compress import Zip

data = 'MULTIPOLYGON (((40 40, 20 45, 45 30, 40 40)), ((20 35, 10 30, 10 10, 30 5, 45 20, 20 35), (30 20, 20 15, 20 25, 30 20)))'
vector = VectorComposite(
    WKT(data),
    ShapefileCompressed(Shapefile(), compress_engine=Zip())
)
  • 从 wkt 读取所有几何图形
vector.geometry_collection()
  • 遍历所有几何
for geom in vector.geometries():
    print(geom)
  • 以 zip 格式创建 shapefile
vector.write('output.zip')
从几何中搜索 UTM 区域
  • 此功能将从某些几何图形中搜索 UTM 区域。
from vectorio.vector import UTMZone, VectorReprojected, WKT
data = 'POLYGON((-73.79131452179155 -11.78691590735885,-27.12139264679149 -12.645910804419744,-47.46330883419978 10.894322081983276,-73.79131452179155 -11.78691590735885))'
ds_wkt = VectorReprojected(WKT(data), out_srid=4326).datasource()
utm = UTMZone()
utm.zones(ds_wkt)  # getting all UTM Zones that intersect with the geometry
  • zone_from_biggest_geom方法只获取一个具有最大几何形状的区域。例如,如果大型几何图形位于 UTM Zone 25N 或 26N,此方法将计算几何图形的面积(对于多边形)或长度(对于线),并返回面积或长度最大的 UTM Zone。

  • 签名zone_from_biggest_geom(self,inp_ds:数据源,wkt_prj_for_metrics:str,in_wkt_prj=PRJ_WGS84)

    • wkt_prj_for_metrics:必需
    • inp_ds : 必需
    • in_wkt_prj:可选。用例输入数据源没有 CRS。
# prj_for_metrics - necessary for make metrics calculations
prj_for_metrics = 'PROJCS["Brazil / Albers Equal Area Conic (WGS84)",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["longitude_of_center",-50.0],PARAMETER["standard_parallel_1",10.0],PARAMETER["standard_parallel_2",-40.0],PARAMETER["latitude_of_center",-25.0],UNIT["Meter",1.0]]'
utm.zone_from_biggest_geom(ds_wkt, wkt_prj_for_metrics=prj_for_metrics) # getting one UTM Zone
  • 在此计算中忽略点。对于仅由点组成的几何图形(多点、点 e 的几何集合等),不建议使用此方法,因为返回的 UTM 区域基于面积/长度度量。很快,将实现一种仅用于点几何的获取 UTM 区域的方法。
具有拓扑错误的几何图形的 UTM 区域

如果您的几何图形存在拓扑错误,您应该使用方法UTMZone().zones_from_iteration(wkt)将 WKT 作为参数传递。此方法将返回 wkt 相交的一组 utm 区域。

utm.zones_from_iteration(
    'POLYGON((-59.408117902654105 -6.5855114909234445,-62.132727277654105 -8.241414689294965,-58.836828840154105 -8.545736525537883,-62.835852277654105 -5.230474986041972,-58.836828840154105 -6.4981931392341705,-59.408117902654105 -6.5855114909234445))'
)
# {'20SW', '21SW'}

引发 Gdal 警告的例外情况

要使用 gdal 警告中的异常,应使用 vectorio.gdal包上的装饰器gdal_warning_as_exception呈现。当使用geometry()方法中的IsValid()方法时,此装饰器将抛出错误。

from vectorio.gdal import gdal_warning_as_exception
from vectorio.vector import WKT

self_intersect_polygon = 'POLYGON((-54.24438490181399 -5.466896872158672,-54.84863294868899 -5.882330540835073,-54.09057630806399 -5.8714019542356475,-54.83764662056399 -5.379399666352095,-54.24438490181399 -5.466896872158672))'

@gdal_warning_as_exception
def possible_error():
    wkt = WKT(self_intersect_polygon)
    ds = wkt.datasource()
    lyr = ds.GetLayer(0)
    feat = lyr.GetFeature(0)
    feat.geometry().IsValid()

possible_error()
# >>> GDALSelfIntersectionGeometry: Self-intersection at or near point -54.469636435829948 -5.6217621987992636
可能的例外
  • GDALSelfIntersectionGeometry:当多边形包含自相交时抛出异常。
  • GDALBadClosedPolygon:当多边形未正确关闭时引发异常。
  • GDALUnknownException:发生未知错误时抛出异常。

Obs:所有异常都在包vectorio.exceptions上可用

直接重投影数据源

  • 要进行空间重投影,请使用上面提到的类 VectorIOReprojected 。但是,在同一时刻将有必要直接重新投影数据源,为此使用 DataSourceReprojected 类。

  • 签名DataSourceReprojected(inp_ds:DataSource,in_srid:int=None,out_srid:int=None,in_wkt_prj=None,out_wkt_prj=None,use_wkt_prj:bool=False)

    • inp_ds:必需。将被重新投影的数据源。
    • in_srid:可选。用例输入数据源没有 CRS。(仅在标志use_wkt_prjFalse时使用)。
    • out_srid:需要use_wkt_prj是否为False
    • in_wkt_prj:可选。OGC 模式上的 WKT 投影。用例输入数据源没有 CRS。(仅当标志use_wkt_prjTrue时使用)。
    • out_wkt_prj:需要use_wkt_prj是否为True。OGC 模式上的 WKT 投影。
    • use_wkt_prj:布尔值,用于定义数据源是否将使用 WKT 投影使用 SRID 重新投影。
  • 方法:

    • ref() -> DataSource:返回重新投影的 DataSource。
  • 使用 SRID 重新投影。

from vectorio.vector import WKT, DataSourceReprojected
data = 'POLYGON((-45.522540331834634 -6.851736627227062,-47.016680956834634 -7.7670786428296275,-45.434649706834634 -8.332736100352385,-45.522540331834634 -6.851736627227062))'
wkt = WKT(data)
ds = wkt.datasource()
new_ds = DataSourceReprojected(ds, out_srid=31983).ref()
wkt.geometry_collection(ds=new_ds)
  • 使用 WKT 投影重新投影。
from vectorio.vector import WKT, DataSourceReprojected
data = 'POLYGON((-45.522540331834634 -6.851736627227062,-47.016680956834634 -7.7670786428296275,-45.434649706834634 -8.332736100352385,-45.522540331834634 -6.851736627227062))'
wkt = WKT(data)
prj = 'PROJCS["SIRGAS 2000 / UTM zone 23S",GEOGCS["SIRGAS 2000",DATUM["Sistema_de_Referencia_Geocentrico_para_America_del_Sur_2000",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6674"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4674"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],AUTHORITY["EPSG","31983"],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
ds = wkt.datasource()
new_ds = DataSourceReprojected(ds, out_wkt_prj=prj, use_wkt_prj=True).ref()
wkt.geometry_collection(ds=new_ds)

计算顶点

所有向量类都有计算所有顶点的方法。下面,以 WKT 为例,但是,此方法也适用于 Shapefile、GeoJson。

  • 签名vertices_by_feature(ds: DataSource) -> dict

    • ds:必填。DataSource 用于计算您的顶点。
  • 返回

    • 将按特征返回一个顶点字典。
from vectorio.vector import WKT
data = 'POLYGON((-45.522540331834634 -6.851736627227062,-47.016680956834634 -7.7670786428296275,-45.434649706834634 -8.332736100352385,-45.522540331834634 -6.851736627227062))'
wkt = WKT(data)
wkt.vertices_by_feature() # works with shapefile or Geojson 
# shapefile_obj.vertices_by_feature() or geojson_obj.vertices_by_feature()

项目详情