從Python使用Python的矩形

Question

我有一個充滿點的文本文件。它們在每行上由逗號限制（x，y）對分隔。例如。

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

我現在需要圍繞每個點創建一個多邊形。多邊形必須具有15公里的緩衝區。我無法訪問ArcGIS或任何其他爲我提供此功能的GIS，因此此時，我想知道是否有人可以幫助我入門？

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())