Rectángulos de puntos usando Python

Question

Tengo un archivo de texto lleno de puntos. Están separados en cada línea por un par de coma limitado (x, y). p.ej.

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

Ahora necesito crear un polígono alrededor de cada uno de estos puntos. El polígono tiene que tener un buffer de 15 kilómetros desde el punto. No tengo acceso a ArcGIS ni a ningún otro SIG que me brinde esta funcionalidad, así que, en este momento, me pregunto si alguien tiene los cálculos que me ayudarán a comenzar.

Answer 1

Quiere usar GDAL/OGR/OSR, que hace proyecciones, búferes, e incluso podría escribir el Shapefile para usted.

Para convertir grados lat/long en metros para su búfer métrico, necesita un sistema de coordenadas proyectado. En el ejemplo siguiente, uso zonas UTM, que se cargan y almacenan en caché dinámicamente. Esto calculará 15 km en el geoide.

También calculo tanto el búfer GIS, que es un polígono redondeado, como el sobre del búfer calculado, que es el rectángulo que busca.

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