I have sets of lat/lon points for vortices in the ocean, and what I want to do is draw curves around them (they are roughly elliptical), and determine their eccentricity and orientation as they move through the ocean in time. I'm using python, and I'm on a mac.
I am having a really dreadful time navigating the sea (pun intended) of special libraries that are intended for this purpose.
Can anyone provide me with a simple list of instructions to do this?
I found some tools that point to Fiona and Shapely as easy-to-use tools, but I cannot for the life of me figure out how to get this array of lat-lon points into a .shp file which Fiona can open.
1) How do I create a .shp file out of these data points? If there is another filetype I should use, please let me know how. This is as far as I've gotten:
import Shapefile
w = shapefile.Writer(shapeType = 8) #8 is multipoint
From here I'm totally stuck, and can't find documentation that helps me past this. I don't need any additional features, literally just the points.
I have tried the following as well: w = shapefile.Writer()
for i in range(vortex.shape[0]):
w.point(vortex[i, 0], vortex[i, 1])
for i in range(0, 327):
print w.shapes()[i].points
w.save('sample_vortex.shp')
Vortex is a numpy array with dimensions (327, 2), with each row containing a lat and a lon value (duh).
This prints out a list of the correct points, so I'm on the right track. I found this snippet of code, which I think should plot my points:
import fiona
import shapely.geometry as geometry
input_shapefile = 'sample_vortex.shp'
shapefile = fiona.open(input_shapefile)
points = [geometry.shape(point['geometry'])
for point in shapefile]
import pylab as pl
x = [p.coords.xy[0] for p in points]
y = [p.coords.xy[1] for p in points]
pl.figure(figsize=(10,10))
pl.plot(x,y,'o', color='#f16824')
But it actually does not, and yields a blank plot, which tell me, I think, that my .shp file isn't correct.
Any ideas?
2) Do you know how to draw these ellipses and get the eccentricity (and how much the major axis deviates from horizontal, while we're at it)?
dlopen(/Users/ben/anaconda/lib/python2.7/site-packages/osgeo/_gdal.so, 2): Library not loaded: @rpath/libjpeg.8.dylib Referenced from: /Users/ben/anaconda/lib/libgdal.20.dylib Reason: image not found