标签: shapely

我有 WGS84 中的 lat、lng 修复列表，我想在这些列表上进行计算，例如点、多边形之间的距离测量等...为此，我计划使用匀称但后来，我需要将其转换为笛卡尔空间，这只是局部准确。

我的问题是我的位置修复可能来自世界各地，所以如果我使用针对我所在地区优化的固定投影，我会在世界其他地方引入错误。是否可以根据当前位置列表的平均位置来定义我自己的以位置对为中心的笛卡尔投影？我需要计算的位置修复总是彼此接近，但不同的位置修复列表可以分布在世界各地。

例如：假设我得到 5 个修复，我需要对其进行计算。然后，我想定义一个在这些定位点附近准确的投影，因为这些定位点之间的距离始终在几公里之内。当我获得接下来的 5 个修复时，它们可能位于世界的完全不同的地方，我想定义一个针对这些位置修复优化的投影。

我将如何解决这个问题？似乎使用 pyproj（如果我理解得很好，则使用 proj.4）是一个好主意，但我无法理解初始化投影所需的字符串，如下所示。有人能帮我吗？

local_proj = pyproj.Proj(r'+proj=tmerc +lat_0=51.178425 +lon_0=3.561298 +ellps=GRS80 +units=meters')

我有以下问题。我有一个装满坐标和三个点的盒子，它们构成了一条线。现在我想计算所有框坐标到该线的最短距离。我有三种方法可以做到这一点，vtk 和 numpy 版本总是有相同的结果，但不是 shapely 的距离方法。但我需要匀称的版本，因为我想测量从一个点到整条线而不是到单独的线段的最近距离。这是迄今为止的示例代码。所以问题是“pdist”：

from shapely.geometry import LineString, Point
import vtk, numpy as np
import itertools

import math

from numpy.linalg import norm

x1=np.arange(4,21)
y1=np.arange(4,21)
z1=np.arange(-7,6)

linepoints = np.array([[1,10,0],[10,10,0],[15,15,0]])


for i in itertools.product(x1,y1,z1):
    
    for m in range(len(linepoints)-1):
        
        line3 = LineString([linepoints[m],linepoints[m+1]])
        
        p = Point(i)
        
        d = norm(np.cross(linepoints[m]-linepoints[m+1], linepoints[m]-i))/norm(linepoints[m+1]-linepoints[m])
        
        dist=math.sqrt(vtk.vtkLine().DistanceToLine(i,linepoints[m],linepoints[m+1]))
        
        pdist = p.distance(line3)
        
        print(d,dist,pdist)

所以我有这样一种情况，我有很多断路的线串，我需要使用 Shapely 的 LineMerge 或 Union OR PostGIS ST_Union 将它们联合在一起。

我现在的想法是使用 Shapely 将Linestrings导入为 Geometry 类型。使用 Shapely 将它们联合或合并，然后导出回数据库中的结果表。

但是，PostGIS 数据库中的几何类型只是一堆乱码。喜欢...

01020000020e61000....

如何使用 Shapely 将其从数据库转换为 Python 几何类型，进行一些操作，然后将其导出回数据库？

目前这是我的代码，它现在只是从数据库中导入 geom 对象字符串并抛出错误，因为它不是 Geometry 类型。

def create_shortest_route_geom(shortest_routes):
    conn = connect_to_database()
    cur = conn.cursor()
    shortest_route_geoms = []
    for route in shortest_routes:
        source = str(int(route[1]))
        target = str(int(route[2]))
        query = 'SELECT the_geom FROM public.ways WHERE target_osm = ' + target + ' AND source_osm = ' + source + ' OR target_osm = ' + …

我可以通过以下方式创建多边形：

#!/usr/bin/env python

from shapely.geometry import Polygon

area = Polygon(((52, 13), (57, 14), (58, 12)))

with open('test.svg', 'w') as f:
    f.write(area.svg())

返回

<path fill-rule="evenodd" fill="#66cc99" stroke="#555555" stroke-width="2.0" opacity="0.6" d="M 52.0,13.0 L 57.0,14.0 L 58.0,12.0 L 52.0,13.0 z" />

这不是有效的 SVG 文件。如何获得有效的 SVG？

我试过

#!/usr/bin/env python

from shapely.geometry import Polygon

area = Polygon(((52, 13), (57, 14), (58, 12)))

with open('test.svg', 'w') as f:
    f.write('<svg version="1.1" xmlns="http://www.w3.org/2000/svg" xmlns:xlink= "http://www.w3.org/1999/xlink">')
    f.write(area.svg())
    f.write('</svg>')

当我查看这个时，视口对于多边形来说太大了。使用 Inkscape 手动编辑它并调整它的大小可以得到：

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<svg
   xmlns:dc="http://purl.org/dc/elements/1.1/"
   xmlns:cc="http://creativecommons.org/ns#"
   xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#" …

我有一个功能

get_polygon(polygon_collection, point):
    for polygon in polygon_collection:
        if polygon.intersects(point):
            return polygon
    return None

这种方法有效，但它在 O(n) * O（单多边形检查）中。如果构建树数据结构，这肯定可以减少到 O(log(n)) * O（单多边形检查）。

匀称是否直接支持？

多边形集合的结构

• 该集合包含 n 个 MultiPolygons - 因此一个对象可以由许多多边形组成。
• 每个 MultiPolygon 可以有飞地/飞地
• 通常 (> 95%)，它们是没有孔的简单多边形
• 通常（> 50％），形状相当简单（例如正方形）

用例示例

多边形列表可以是德国的邮政编码区域。那将是几千。然后我有我和一些朋友的 GPS 位置，也有几千个。我想说我们在哪个区域获得了最多的数据点。

在对多边形进行分割等操作之后，我想验证它是否是一个矩形。
我试过simplify然后数数是否coords是5...

>>> from shapely.geometry import Polygon
>>> from shapely.ops import split
>>> 
>>> poly1 = Polygon([(0, 0), (0, 1), (0, 3), (2, 3), (2, 2), (2, 0), (0, 0)])
>>> 
>>> poly_check=poly1.simplify(0)
>>> if len(poly_check.exterior.coords)==5:
>>>     print 'Yes, it is a rectangle...'
>>> else:
>>>     print 'No, it is not a rectangle...'
>>> 
Yes, it is a rectangle...

但如果起点位于边缘的中间，则该方法不起作用。

>>> #poly2 is actually a rectangle
>>> poly2 = Polygon([(0, 1), (0, 3), (2, 3), (2, …

我有一个简单的问题，但我找不到答案我正在寻找“最小旋转矩形”多边形相对于纬度或经度的主轴角

df4.minimum_rotated_rectangle

有人有这个库存吗 提前致谢

我正在开发一个需要坐标映射的项目 - 确定坐标点是否存在于一系列多边形中。映射的数量相当大 - 跨越 100 多个多边形的约 1000 万个坐标。

在继续之前，我已经查看了此处和此处的问题。这个问题并不多余，因为它涉及动态点和静态多边形。

我通过在 200 万个多边形的子集中映射单个坐标来缩小该问题的项目范围。这是我使用的代码：

from shapely.geometry import shape, Point

f = open('path/to/file.geojson', 'r')
data = json.loads(f.read())

point = Point(42.3847, -71.127411)
for feature in data['features']:
    polygon = shape(feature['geometry'])
    if polygon.contains(point):
        print(polygon)

迭代 200 万个多边形（在本例中为建筑足迹）大约需要 30 秒（时间太长）。

我也尝试过使用mplPath如下：

import matplotlib.path as mplPath

building_arrays = [np.array(data['features'][i]['geometry']['coordinates'][0])
                   for i, v in enumerate(tqdm(data['features']))]
bbPath_list = [mplPath.Path(building)
               for building in tqdm(building_arrays)]

for b in tqdm(bbPath_list):
    if b.contains_point((-71.1273842, 42.3847423)):
        print(b)

这大约需要 6 秒。一个改进，但考虑到我需要的映射量，仍然有点慢。 …

我正在使用此代码，并且在另一个窗口中打开了一个可视化/图表，在 ipython shell 中没有任何问题：

In [1]: import numpy as np                                                                                                                            

In [2]: import matplotlib.pyplot as plt                                                                                                               

In [3]: matplotlib                                                                                                                                    
Using matplotlib backend: Qt5Agg

In [4]: x = np.linspace(0, 3*np.pi, 500)                                                                                                              

In [5]: plt.plot(x, np.sin(x**2))                                                                                                                     
Out[5]: [<matplotlib.lines.Line2D at 0x7fb6d69ab470>]

In [6]:  

检查此屏幕截图。如果我使用 matplotlib 以外的任何库，那么我不会得到任何可视化/图表。我确实在笔记本中得到了它，但在 ipython shell 中却没有。难道我做错了什么 ？

In [7]: poly = Polygon([(0,0), (0,5), (5,5), (5,0)])                                                                                                  

In [8]: print(poly)                                                                                                                                   
POLYGON ((0 0, 0 5, 5 5, 5 0, 0 0))

In [9]: print('area', poly.area)                                                                                                                      
area 25.0

In [10]: display(poly) …

我有两个一维数组，想将它们组合成一个 Point GeoSeries，如下所示：

import numpy as np
from geopandas import GeoSeries
from shapely.geometry import Point

x = np.random.rand(int(1e6))
y = np.random.rand(int(1e6))
GeoSeries(map(Point, zip(x, y)))

在我的笔记本电脑上大约需要 5 秒。是否可以加速GeoSeries的生成？