[GIS] Buffer miles in decimal degrees with fiona and shapely

coordinate systemfionapythonshapely

I have a script where I parse geojson earthquake data from around the world and create the point locations of the earthquake incidents. The script works fine. The CRS is 4326 or wgs84/lat long. I am planning on creating certain distances in miles around the points based on their magnitudes. For testing and learning purposes I want to create a 50 mile buffer around each point. I am using fiona and shapely. I can post all my code if somebody asks, although it does not seem necessary

import fiona
import shapely
python code.....
crs = from_epsg(4326)
point = Point(float(i["geometry"]["coordinates"][0]),float(i["geometry"]["coordinates"][1]))
output.write({'properties':{'Place': place, 'Magnitude': mag, 'KM': km, 'Bearing': bear},'geometry': mapping(shape(point.buffer(5)))})

I have seen this post
Create a buffer in decimal degrees

but I do not want to call in arcpy for this analysis.

I am open to a custom Python conversion from miles to decimal degrees or even feet to decimal degrees.

Best Answer

You can create a custom transverse mercator projection on the earthquake center, draw a 50 miles circle in that CRS, and reproject it to EPSG:4326.

This is handsome for automation, because you can always use the same circle for all earthquakes.

Related Solutions

[GIS] problems with py2exe, shapely and fiona

I had the same issue when packaging a fiona project with py2exe. I think py2exe simply cannot resolve all the fiona dependencies. But fortunately fiona is well behaved: just add it under packages into the py2exe script and it will probably also work for you:

opts = {
  'py2exe': {
    'packages':['fiona'],
    'bundle_files': 3
    }
  }

setup(
      console=[{"script" : "test.py"}],
      options=opts,
      zipfile=None
      )

[GIS] Find common borders of features and get their attributes using fiona and shapely

As Fiona works with dictionaries, it is easy.

With the example of Finding the common borders between polygons in the same shapefile

import fiona
from shapely.geometry import shape, mapping
# load the shapefile as an iterable
polys = fiona.open("polygons.shp")

Now, with the itertools module:

import itertools
for pol1,pol2 in  itertools.combinations(polys, 2):
    if shape(pol1['geometry']).touches(shape(pol2['geometry'])):
        print  "Polygon ",pol1['properties']['id'],  " touches Polygon " , pol2['properties']['id']

Polygon 1 touches Polygon 4
Polygon 2 touches Polygon 4
Polygon 3 touches Polygon 4

And directly:

schema = {'geometry': 'LineString','properties':  OrderedDict([(u'poly1', 'int:10'), (u'poly2', 'int:10')])}
with fiona.open('intersection2.shp','w','ESRI Shapefile', schema) as output:
   with fiona.open("intersections.shp") as polys:
      for pol1,pol2 in  itertools.combinations(polys, 2):
          geompol1 = shape(pol1['geometry'])
          geompol2 = shape(pol2['geometry'])
          if geompol1.touches(geompol2):
             output.write({'geometry':mapping(geompol1.intersection(geompol2)),  'properties': {'poly1':pol1['properties']['id'],'poly2':pol2['properties']['id']}})

Control

 for line in fiona.open("intersection2.shp"):
     print line['properties']
 OrderedDict([(u'poly1', 1), (u'poly2', 4)])
 OrderedDict([(u'poly1', 2), (u'poly2', 4)])
 OrderedDict([(u'poly1', 3), (u'poly2', 4)])

Related Question