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.

Community
  • 1
ziggy
  • 4,515
  • 3
  • 37
  • 83
  • 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 00:22
  • What do you mean by geometry service – ziggy Nov 16 '16 at 00:25
  • So is it even possible to apply 50 miles buff zones to these points? – ziggy Nov 16 '16 at 00:26
  • A service that takes geometry requests and returns JSON coordinate arrays. – Vince Nov 16 '16 at 00:27
  • 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 00:29
  • 1
    Reproject points to appropriate CRS, buffer, then reproject buffer polys back to EPSG:4326 if desired - http://gis.stackexchange.com/a/127432/2856 – user2856 Nov 16 '16 at 00:31

1 Answers1

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.

AndreJ
  • 76,698
  • 5
  • 86
  • 162
  • 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 http://gis.stackexchange.com/questions/83861/using-customized-coordinate-system-for-archaeological-site-data, 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
  • so from reading your answer it seems I would have to use the qGIS create custom projection, save that file and call it from my python script? is there a way to define the custom projection in the python script and use it with shapely/fiona? – ziggy Nov 16 '16 at 17:06
  • i know in arcgis you can do something like this: coordinateSystem = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]" – ziggy Nov 16 '16 at 17:08
  • 1
    I don't have experience with python and fiona, but this page tells you how to define a CRS without EPSG code: http://toblerity.org/fiona/manual.html#format-drivers-crs-bounds-and-schema using from_string. – AndreJ Nov 16 '16 at 17:57
  • andre i am having a bit of trouble understanding how to construct the custom CS. when you say "create a custom transverse mercator projection on the earthquake center" and in my custom CS I have to fill in the lat_0 and lon_0 for the center point. 1. do you mean for each earthquake the lat_0 and lon_0 values have to change? 2. this string of a cs {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'} what is no_defs? and how would I go about writing your custom recommendation to this format? – ziggy Nov 19 '16 at 17:12
  • ps: if this is too much of a hassle for you dont worry about it, Im just a bit confused – ziggy Nov 19 '16 at 17:12
  • 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
  • okay i realized it would be easier for the CRS to be in pyproj - what would that syntax look like? also if there is a resource where I can learn about all what the projection detail syntax stands for and means that would be extremely helpful to learn on my own – ziggy Nov 19 '16 at 20:25
  • 1
    The syntax is almost the same: http://gis.stackexchange.com/questions/78838/how-to-convert-projected-coordinates-to-lat-lon-using-python/189695#189695 and http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/. You might look into https://github.com/OSGeo/proj.4/wiki for syntax details. – AndreJ Nov 20 '16 at 07:09