Python Rectangles with Dots

Question

I have a text file full of dots. They are separated on each line using a limited comma (x, y) pair. eg.

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

Now I need to create a polygon around each of these points. The polygon must have a 15 km buffer from the point. I don’t have access to ArcGIS or any other GIS that provides me this function, so at the moment I’m wondering if anyone has any math to help me get started?

Answer 1

You want to use GDAL / OGR / OSR , which executes forecasts, buffers, and can even write a Shapefile for you.

To convert lat / long degrees to meters for your metric buffer, you need a projected coordinate system. In the example below, I use UTM zones that are dynamically loaded and cached. It will calculate 15 km on a geoid.

I also compute both the GIS buffer, which is a rounded polygon, and the envelope from the computed buffer, which is the rectangle you are looking for.

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