Rettangoli da punti usando Python

Question

Ho un file di testo pieno di punti. Sono separati su ogni riga da una coppia limitata (x, y). per esempio.

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

Ora ho bisogno di creare un poligono attorno a ciascuno di questi punti. Il poligono deve avere un buffer di 15 km dal punto. Non ho accesso ad ArcGIS o ad altri GIS che forniscono questa funzionalità per me, quindi a questo punto, mi chiedo se qualcuno ha la matematica che mi aiuterà a iniziare?

Answer 1

Si desidera utilizzare GDAL/OGR/OSR, che esegue proiezioni, buffer e potrebbe persino scrivere lo Shapefile per conto dell'utente.

Per convertire i gradi lat/long in metri per il buffer metrico, è necessario un sistema di coordinate proiettato. Nell'esempio seguente utilizzo le zone UTM, che vengono caricate e memorizzate in modo dinamico. Questo calcolerà 15 km sul geoide.

Inoltre, calcolo sia il buffer GIS, che è un poligono arrotondato, sia la busta dal buffer calcolato, che è il rettangolo che si cerca.

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