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.

  • 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 Nov 16 '16 at 0:22
  • What do you mean by geometry service – ziggy Nov 16 '16 at 0:25
  • So is it even possible to apply 50 miles buff zones to these points? – ziggy Nov 16 '16 at 0:26
  • 1
    Simplest/most accurate is to project point to UTM or similiar CRS, buffer point, project buffer back to 4326. – klewis Nov 16 '16 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 Nov 16 '16 at 0:31
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.

  • 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 Nov 16 '16 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 Nov 16 '16 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 Nov 16 '16 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 Nov 19 '16 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 Nov 20 '16 at 7:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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