Hình chữ nhật từ các điểm sử dụng Python

Question

Tôi có một tệp văn bản đầy đủ các điểm. Chúng được phân cách trên mỗi dòng bằng cặp dấu phẩy (x, y). ví dụ.

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

Bây giờ tôi cần tạo đa giác xung quanh mỗi điểm này. Đa giác phải có bộ đệm 15 km từ điểm. Tôi không có quyền truy cập vào ArcGIS hoặc bất kỳ GIS nào khác cung cấp chức năng này cho tôi nên tại thời điểm này, tôi tự hỏi liệu có ai có toán học sẽ giúp tôi bắt đầu không?

Answer 1

Bạn muốn sử dụng GDAL/OGR/OSR, dự đoán, bộ đệm và thậm chí có thể viết Shapefile cho bạn.

Để chuyển đổi độ vĩ độ/độ dài thành mét cho bộ đệm chỉ số của bạn, bạn cần có hệ tọa độ được chiếu. Trong ví dụ dưới đây tôi sử dụng các vùng UTM, được nạp động và lưu trữ. Điều này sẽ tính toán 15 km trên geoid.

Ngoài ra tôi tính cả bộ đệm GIS, là một đa giác tròn, và phong bì từ bộ đệm được tính toán, là hình chữ nhật mà bạn tìm kiếm.

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