6

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.

Improve this question
6
  • The aren't any fixed conversions between angular and linear units. In the ArcGIS platform, a geometry service would be used to perform the calculation (not arcpy).
    – Vince
    Commented Nov 16, 2016 at 0:22
  • What do you mean by geometry service
    – ziggy
    Commented Nov 16, 2016 at 0:25
  • So is it even possible to apply 50 miles buff zones to these points?
    – ziggy
    Commented Nov 16, 2016 at 0:26
  • 1
    Simplest/most accurate is to project point to UTM or similiar CRS, buffer point, project buffer back to 4326.
    – klewis
    Commented Nov 16, 2016 at 0:29
  • 1
    Reproject points to appropriate CRS, buffer, then reproject buffer polys back to EPSG:4326 if desired - gis.stackexchange.com/a/127432/2856
    – user2856
    Commented Nov 16, 2016 at 0:31

1 Answer 1

Reset to default
3

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.

Improve this answer
10
  • this sounds like a great idea, is there any good links or tutorials you can provide in creating a custom transverse mercator projection- i have never done such a task
    – ziggy
    Commented Nov 16, 2016 at 15:28
  • See gis.stackexchange.com/questions/83861/…, first part of my answer there. You have to fill in the lat_0 and lon_0 for the center point.
    – AndreJ
    Commented Nov 16, 2016 at 16:20
  • 1
    I don't have experience with python and fiona, but this page tells you how to define a CRS without EPSG code: toblerity.org/fiona/… using from_string.
    – AndreJ
    Commented Nov 16, 2016 at 17:57
  • 1
    No poroblem: assuming the earthquake is at 51.4N 7.0E, just use from_string("+proj=tmerc +lat_0=51.4 +lon_0=7 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"). For every other earthquake, exchange the lat_0 and lon_0. You might use a string variable in which you can exchange the values.
    – AndreJ
    Commented Nov 19, 2016 at 19:23
  • 1
    The syntax is almost the same: gis.stackexchange.com/questions/78838/… and all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj. You might look into github.com/OSGeo/proj.4/wiki for syntax details.
    – AndreJ
    Commented Nov 20, 2016 at 7:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.