A pure python implementation of the elliptical Fourier analysis method described by Kuhl and Giardina (1982). This package is designed to allow the rapid analysis of spatial data stored as ESRI shapefiles, handling all of the geometric conversions. The resulting data can be written back to shapefiles to allow analysis with other spatial data or can be plotted using matplotlib.
The code is built upon the pyefd module and it is hoped that this package will allow more geoscientists to apply this technique to analyze spatial data using the elliptical Fourier descriptor technique as there is no longer a data conversion barrier to entry. This package is also more feature rich than previous implementations, providing calculations of Fourier power and spatial averaging of collections of ellipses.
Examples of Fourier ellipses (black) being fitted to a shapefile outline (red), for increasing numbers of harmonics.
- Built-in geometry processing, just pass in a shapefile and get results quickly!
- Fourier coefficient average and standard deviation calculation
- Handles spatial input data through the pyshp library
- Compute an appropriate number of harmonics for a given polygon
- Basic plotting for analysis and debugging through matplotlib
- Write Fourier ellipses as shapefiles
Install spatial_efd by running:
$ pip install spatial_efdThis package supports Python 2.7 and Python 3 and is tested on Linux and Windows environments, using both the standard python interpreter and pypy. It requires matplotlib, numpy, future and pyshp. These packages will all install automatically if spatial_efd is installed using pip.
Note: pypy is currently not being tested via CI due to a matplotlib build error. Please get in touch if this is an issue.
Dependencies can be tracked by visiting requires.io
Note that Python 2 has reached end of life and although the code currently works under Python 2, this will not be supported, and future updates may completely break Python 2 support without warning.
A range of unit tests are included in the /test/ directory. These can be run using pytest:
$ pytestMany of these tests make use of the example_data.shp file which is a shapefile containing six polygons taken from a real dataset of landslide source areas.
The first step in using spatial_efd is always to load a shapefile:
import spatial_efd
shp = spatial_efd.LoadGeometries('test/fixtures/example_data.shp')This creates a shapefile object shp which contains the polygon geometries we want to analyze. As in most cases more than one polygon will be stored in an individual file, a single polygon can be selected for processing using python's list notation:
x, y, centroid = spatial_efd.ProcessGeometryNorm(shp[1])This loads the geometry from the 2nd polygon within the shapefile into a list of x and a list of y coordinates. This method also computes the centroid of the polygon, which can be useful for later analysis. To make comparisons between data from different locations simpler, these data are normalized.
If you already know how many harmonics you wish to compute this can be specified during the calculation of the Fourier coefficients:
harmonic = 20
coeffs = spatial_efd.CalculateEFD(x, y, harmonic)However, if you need to quantify the number of harmonics needed to exceed a threshold Fourier power, this functionality is available. To do this, an initial set of coefficients need to be computed to the number of harmonics required to equal the Nyquist frequency:
nyquist = spatial_efd.Nyquist(x)
tmpcoeffs = spatial_efd.CalculateEFD(x, y, nyquist)
harmonic = spatial_efd.FourierPower(tmpcoeffs, x)
coeffs = spatial_efd.CalculateEFD(x, y, harmonic)Once the coefficients have been calculated they can be normalized following the steps outlined by Kuhl and Giardina (1982):
coeffs, rotation = spatial_efd.normalize_efd(coeffs, size_invariant=True)size_invariant should be set to True (the default value) in most cases to normalize the coefficient values, allowing comparison between polygons of differing sizes. Set size_invariant to False if it is required to plot the Fourier ellipses alongside the input shapefiles, or if the Fourier ellipses are to be written to a shapefile. These techniques which apply to normalized data are outlined later in this document.
A set of coefficients can be converted back into a series of x and y coordinates by performing an inverse transform, where the harmonic value passed in will be the harmonic reconstructed:
xt, yt = spatial_efd.inverse_transform(coeffs, harmonic=harmonic)Wrappers around some of the basic matplotlib functionality is provided to speed up the visualization of results:
ax = spatial_efd.InitPlot()
spatial_efd.PlotEllipse(ax, xt, yt, color='k', width=1.)
spatial_efd.SavePlot(ax, harmonic, '/plots/myfigure', 'png')This example generates an axis object, plots our transformed coordinates onto it with a line width of 1 and a line color of black. These axes are saved with a title denoting the harmonic used to generate the coordinates and are saved in the format provided in the location provided.
Note that as this plotting is performed using matplotlib many other formatting options can be applied to the created axis object, to easily create publication ready plots.
To plot an overlay of a Fourier ellipse and the original shapefile data, a convenience function has been provided to streamline the coordinate processing required.
Plotting the normalized coefficients, where the data has been processed using the ProcessGeometryNorm method is undertaken as follows (Note that size_invariant has been set to False):
# size_invariant must be set to false if a normalized Fourier ellipse
# is to be plotted alongside the shapefile data
coeffs, rotation = spatial_efd.normalize_efd(coeffs, size_invariant=False)
ax = spatial_efd.InitPlot()
spatial_efd.plotComparison(ax, coeffs, harmonic, x, y, rotation=rotation)
spatial_efd.SavePlot(ax, harmonic, '/plots/myComparison', 'png')Which produces a figure like this:
Example of a normalized Fourier ellipse (black) being plotted on top of a shapefile outline (red).
All of the above examples have focused on processing a single polygon from a multipart shapefile, but in most cases multiple geometries will be required to be processed. One of the common techniques surrounding elliptical Fourier analysis is the averaging of a collection of polygons. This can be achieved as follows:
shp = spatial_efd.LoadGeometries('test/fixtures/example_data.shp')
coeffsList = []
for shape in shp:
x, y, centroid = spatial_efd.ProcessGeometryNorm(shape)
harmonic = 10
coeffs = spatial_efd.CalculateEFD(x, y, harmonic)
coeffs, rotation = spatial_efd.normalize_efd(coeffs, size_invariant=True)
coeffsList.append(coeffs)
avgcoeffs = spatial_efd.AverageCoefficients(coeffsList)Once the average coefficients for a collection of polygons has been computed, the standard deviation can also be calculated:
SDcoeffs = spatial_efd.AverageSD(coeffsList, avgcoeffs)With the average and standard deviation coefficients calculated, the average shape, with error ellipses can be plotted in the same manner as individual ellipses were plotted earlier
x_avg, y_avg = spatial_efd.inverse_transform(avgcoeffs, harmonic=harmonic)
x_sd, y_sd = spatial_efd.inverse_transform(SDcoeffs, harmonic=harmonic)
ax = spatial_efd.InitPlot()
spatial_efd.PlotEllipse(ax, x_avg, y_avg, color='b', width=2.)
# Plot avg +/- 1 SD error ellipses
spatial_efd.PlotEllipse(ax, x_avg + x_sd, y_avg + y_sd, color='k', width=1.)
spatial_efd.PlotEllipse(ax, x_avg - x_sd, y_avg - y_sd, color='k', width=1.)
spatial_efd.SavePlot(ax, harmonic, '/plots/average', 'png')Which produces a figure like this:
Example of an average Fourier ellipse (blue) being plotted with standard deviation error ellipses (black).
In cases where the original coordinates are needed, a different processing method can be called when loading coordinates from a shapefile, to return the non-normalized data:
x, y, centroid = spatial_efd.ProcessGeometry(shp[1])This method should be used where the original coordinates need to be preserved, for example if output to a shapefile is desired. To plot non-normalized data alongside the original shapefile data, the locus of the coefficients must be computed and passed as an argument to the inverse transform method:
locus = spatial_efd.calculate_dc_coefficients(x, y)
xt, yt = spatial_efd.inverse_transform(coeffs, harmonic=harmonic, locus=locus)To plot non-normalized coefficients, again call the plotComparison method, with the rotation value set to 0 as no normalization has been performed on the input data:
ax = spatial_efd.InitPlot()
spatial_efd.plotComparison(ax, coeffs, harmonic, x, y, rotation=0.)
spatial_efd.SavePlot(ax, harmonic, '/plots/myComparison', 'png')Which produces a figure like this:
Example of a non-normalized Fourier ellipse (black) being plotted on top of a shapefile outline (red).
In the case of the non-normalized data plotted above, these ellipses can also be written to a shapefile to allow further analysis in a GIS package:
shape_id = 1
shpinstance = spatial_efd.generateShapefile('mydata/myShapefile', prj='test/fixtures/example_data.prj')
shpinstance = spatial_efd.writeGeometry(coeffs, x, y, harmonic, shpinstance, shape_id)The first method called creates a blank shapefile in the path mydata, ready to be populated with Fourier ellipses. By passing in the existing example.prj file to the save method, a new projection file will be generated for the saved shapefile, ensuring that it has the correct spatial reference information for when it is loaded into a GIS package. Note that no reprojection is performed as the aim is for the input and output coordinate systems to match. If this parameter is excluded, the output shapefile will have no defined spatial reference system.
The second method can be wrapped in a loop to write as many ellipses as required to a single file. shape_id is written into the attribute table of the output shapefile and can be set to any integer as a means of identifying the Fourier ellipses.
For more detailed guidance on all of the functions and arguments in this package please check out the source code on github or the API documentation.
I welcome contributions to the code, head to the issue tracker on GitHub to get involved!
If you find any bugs, have any questions or would like to see a feature in a new version, drop me a line:
- Twitter: @GIStuart
- Email: stuart@swdg.io
The project is licensed under the MIT license.
If you use this package for scientific research please cite it as:
Grieve, S. W. D. (2017), spatial-efd: A spatial-aware implementation of elliptical Fourier analysis, The Journal of Open Source Software, 2 (11), doi:10.21105/joss.00189.
You can grab a bibtex file here.
Kuhl and Giardina (1982). Elliptic Fourier features of a closed contour. Computer graphics and image processing, 18(3), 236-258.