问题描述
我有一个充满分数的文本文件.它们在每行上用逗号限制(x,y)对分隔.例如.
I have a text file full of points. They are separated on each line by a comma-limited (x,y) pair. eg.
-43.1234,40.1234\n
-43.1244,40.1244\n
etc.
我现在需要围绕每个点创建一个多边形.从该点开始,多边形必须有15公里的缓冲区.我无权访问为我提供此功能的ArcGIS或任何其他GIS,因此,我现在想知道是否有人掌握了可以帮助我入门的数学方法?
I now need to create a polygon around each of these points. The polygon has to have a 15 kilometer buffer from the point. I don't have access to ArcGIS or any other GIS that provides this functionality for me so at this point, I am wondering if anyone has the math that will help me get started?
推荐答案
您要使用 GDAL/OGR/OSR ,它可以进行投影,缓冲,甚至可以为您编写Shapefile.
You want to use GDAL/OGR/OSR, which does projections, buffers, and could even write the Shapefile for you.
要为公制缓冲区将纬度/经度转换为米,您需要一个投影坐标系.在下面的示例中,我使用动态加载和缓存的UTM区域.这将在大地水准面上计算出15公里.
In order to convert degrees lat/long into metres for your metric buffer, you need a projected coordinate system. In the below example I use UTM zones, which are dynamically loaded and cached. This will calculate 15 km on the geoid.
我还要计算GIS缓冲区(一个圆形多边形)和计算缓冲区中的包络(您要寻找的矩形).
Also I calculate both the GIS buffer, which is a roundish polygon, and the envelope from the calculated buffer, which is the rectangle you seek.
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())
这篇关于使用Python从点矩形的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!