Polygon Orientation Property and Area Calculation

[santisoler]

On this PR I divide the Polygon improvements with the ones made in Moulder on #411 .
I've cleaned up the code in a more Pythonic way, adding properties, property setters decorators and private methods.

In summary, this PR adds:

• Area computation of the Polygon.
• Get the Polygon orientations through the area computation without applying the absolute value. If area is positive, the orientation is counterclockwise. If area is negative, the orientation is clockwise.
• Change the Polygon orientation by rearranging its the vertices.
• Two new tests: one for the area computation (comparing the code with a regular hexagon) and one for the orientation.

I've founded out that the Polygon orientation (clockwise or counterclockwise) depends on the direction of the cartesian axes on which it is defined.

For example, in a x, z system (with x->Right and z->Down, as we assume in talwani.gz), a clockwise Polygon has a negative area, so the code interprets it correctly as clockwise.

On the other hand, in a regular cartesian x, y system (with x->Right, y->Up), a clockwise Polygon has positive area, but the code will interpret it as a counterclockwise Polygon.

In the following scripts I intent to show this behaviour:

from __future__ import division
import numpy as np
from fatiando.mesher import Polygon

angles = np.arange(0, 360, 60)
x, z = np.cos(angles*np.pi/180), np.sin(angles*np.pi/180)
vertices = zip(x, z)
poly = Polygon(vertices, force_clockwise=False)
print poly.orientation
print poly._calculate_area(absolute=False)

>> clockwise
>> -2.59807621135

from __future__ import division
import numpy as np
from fatiando.mesher import Polygon

angles = np.arange(0, 360, 60)[::-1]
x, y = np.cos(angles*np.pi/180), np.sin(angles*np.pi/180)
vertices = zip(x, y)
poly = Polygon(vertices, force_clockwise=False)
print poly.orientation
print poly._calculate_area(absolute=False)

>> counterclockwise
>> 2.59807621135

I can conclude that a positive area implies that the Polygon orientation follows the same right-hand rule than the axes. If area is negative, the Polygon orientation is contrary to the right-hand rule applied to the axes (hope you understand this conclusion).

Therefore the clockwise or counterclockwise orientation is defined by the pair of axes we use.
Should we stick to only one type of axes or give the option to choose between a right or left handed axes?
What's the best way to do this?

Checklist:

• Make tests for new code (at least 80% coverage)
• Create/update docstrings
• Include relevant equations and citations in docstrings
• Docstrings follow the style conventions
• Code follows PEP8 style conventions
• Code and docs have been spellchecked
• Include new dependencies in doc/install.rst, environment.yml, ci/requirements.txt.
• Documentation builds properly (run cd doc; make locally)
• Changelog entry (leave for last to avoid conflicts)
• First-time contributor? Add yourself to doc/contributors.rst (leave for last to avoid conflicts)

[leouieda]

@santis19 thanks for raising this issue with the coordinates. I never liked the hacky way we dealt with it (pretty much just ignored it).

The modeling code will need to know the orientation of the polygon (horizontal or vertical). And it doesn't make sense to use "x" and "y" for vertices coordinates in vertical polygons. So I think the best would be to have two classes: PolygonHorizontal and PolygonVertical. They can both inherit from a common private class BasePolygon that implements most of the mechanics in the current Polygon. The two classes would only need to chance the attribute "y" for "z" and set their orientation.

Also remember that x->north and y->east.

[santisoler]

@leouieda, I was thinking the Polygon class without any specific application, just like a geometrical element independent of how it will be used.

If on the other hand we take into account our purposes, the Polygon class is intended to be used in the forward gravity calculations through the Talwani method.
With this in mind, would a PolygonHorizontal class ever be used? Remember that Talwani equations are meant to be used on computations point belonging to the same plane than the polygon (i.e. vertical polygons).

[leouieda]

The only other place I think of a horizontal polygon being useful is for the PolygonalPrisms. They have a topolygon method that was used for plotting. It's kind of silly because anything that can plot a Polygon could probably plot a PolygonalPrism without any changes.

You're right, we can have just the vertical polygon. It could be the default Polygon and if we ever need it we can have a PolygonHorizontal later. You might need to delete the topolygon method and see if it's being used anywhere.

[santisoler]

I've modified the Polygon class in order to assume it as a vertical polygon with clockwise orientation as default.

• I've changed the y to z in the vertices' coordinates.
• I've added a third orientation status: undefined; for area == 0 cases (I don't know if this could happen, but I've added it as a safeguard).
• I've modified the talwani.gz() function for the new Polygon. It will ignore non clockwise oriented Polygons and will raise a warning.
• I've also modified the documentation of Polygon and talwani.
• I've added Polygon tests, specially for area and orientation calculation.
• I've deleted the topolygon method of PolygonalPrism.

Furthermore, I've seen that Sphere is a subclass of Polygon. I think it deserves a few changes, specially in the documentation.

[santisoler]

@leouieda I've noticed that if we adequate Square to work with the same coordinate system as the vertical Polygon, a lot of errors involving the fatiando.seismic module raise.

Maybe just for the Square case we can define a vertical one and an horizontal one. The first one intended to be just a particular case of Polygon, while the other one will be intended to work with the fatiando.seismic module.

I did't get into the fatiando.seismic and don't know if it always uses vertical squares. In such case, a single vertical Square will work, and just a few changes in fatiando.seismic variable names must be carried on.

[leouieda]

@santis19 Square is used in the forward modeling for srtomo (SquareMesh actually). And that kind of assumes that it's horizontal. Though maybe it shouldn't? This raises the point that sometimes we'll need horizontal polygons, if only so that we can have horizontal squares.

I'd say, leave Square as it is right now (being horizontal) and add the PolygonHorizontal class. Open an issue for the Square problem and we'll work on it in a separate PR.

[santisoler]

@leouieda I've added two new classes: PolygonBase and PolygonHorizontal.

• PolygonBase is intended to be the base class for Polygon PolygonHorizontal.

• PolygonHorizontal defines horizontal polygons, and Square class is now a subclass of it.

• Polygon class defines vertical polygons. I didn't rename the class in order to keep compatibility with the current stable release. Is this ok?

With this definition of Square, fatiando.seismic doesn't break. Maybe some changes on its documentation is needed. Can we say there's no issue regarding srtomo?

[leouieda]

@santis19 sounds great! There doesn't seem to be any problem with srtomo.

I'll update your branch so that you stop getting these failures on TravisCI. After that, this is good to merge.