使用Python从点生成矩形

Question

我有一个文本文件，里面存储了很多点的信息。每行一对(x,y)坐标值，用逗号分隔。例如：

-43.1234,40.1234\n
-43.1244,40.1244\n
etc.

我现在需要在每个点周围创建一个多边形。 多边形必须从该点向外缓冲15公里。 我没有访问ArcGIS或任何其他为我提供此功能的GIS，因此此时，我想知道是否有人拥有可帮助我入门的数学知识？

- aeupinhere

---

你是否正在寻找与ArcGIS兼容的文件格式？ - Usagi

此外，您是否处于投影坐标系中？ 您需要处于投影坐标系中，以便使用缓冲区获取准确的结果。 您可能希望尝试在GIS Stack Exchange网站上提出这个问题。 - Usagi

我只是将一切写入到shapefile中。我已经知道如何做了，但这不是很“Pythonic”。http://forums.esri.com/Thread.asp?c=93&f=983&t=289084。我正在寻找更加优美的解决方法。 - aeupinhere

数据的格式是什么？是经度/纬度、纬度/经度、X/Y、Y/Z等吗？ - Mike T

Answer 1

你想使用 GDAL/OGR/OSR，它可以进行投影、缓冲区处理，甚至可以为你编写 Shapefile 文件。

为了将经纬度转换为米数以进行度量缓冲，你需要一个投影坐标系。在下面的示例中，我使用 UTM 带，它们是动态加载和缓存的。这将在大地椭球上计算出 15 公里。

此外，我还计算了 GIS 缓冲区，它是一个圆形多边形，以及从计算缓冲区得出的包络矩形，这是你要寻找的矩形。

from osgeo import ogr, osr

# EPSG:4326 : WGS84 lat/lon : http://spatialreference.org/ref/epsg/4326/
wgs = osr.SpatialReference()
wgs.ImportFromEPSG(4326)    
coord_trans_cache = {}
def utm_zone(lat, lon):
    """Args for osr.SpatialReference.SetUTM(int zone, int north = 1)"""
    return int(round(((float(lon) - 180)%360)/6)), int(lat > 0)

# Your data from a text file, i.e., fp.readlines()
lines = ['-43.1234,40.1234\n', '-43.1244,40.1244\n']
for ft, line in enumerate(lines):
    print("### Feature " + str(ft) + " ###")
    lat, lon = [float(x) for x in line.split(',')]
    # Get projections sorted out for that UTM zone
    cur_utm_zone = utm_zone(lat, lon)
    if cur_utm_zone in coord_trans_cache:
        wgs2utm, utm2wgs = coord_trans_cache[cur_utm_zone]
    else: # define new UTM Zone
        utm = osr.SpatialReference()
        utm.SetUTM(*cur_utm_zone)
        # Define spatial transformations to/from UTM and lat/lon
        wgs2utm = osr.CoordinateTransformation(wgs, utm)
        utm2wgs = osr.CoordinateTransformation(utm, wgs)
        coord_trans_cache[cur_utm_zone] = wgs2utm, utm2wgs
    # Create 2D point
    pt = ogr.Geometry(ogr.wkbPoint)
    pt.SetPoint_2D(0, lon, lat) # X, Y; in that order!
    orig_wkt = pt.ExportToWkt()
    # Project to UTM
    res = pt.Transform(wgs2utm)
    if res != 0:
        print("spatial transform failed with code " + str(res))
    print(orig_wkt + " -> " + pt.ExportToWkt())
    # Compute a 15 km buffer
    buff = pt.Buffer(15000)
    print("Area: " + str(buff.GetArea()/1e6) + " km^2")
    # Transform UTM buffer back to lat/long
    res = buff.Transform(utm2wgs)
    if res != 0:
        print("spatial transform failed with code " + str(res))
    print("Envelope: " + str(buff.GetEnvelope()))
    # print("WKT: " + buff.ExportToWkt())

- Mike T