3

I've download from Mapzen many geoJson files related to countries with their respective administration level. Then I've parsed all the geoJson files and saved them to my sql server 2012 as geography spatial data.

I've then done a check of the Luxemburg map and noticed some strange results. Some of the maps where created correctly, with the corresponding inner area:

enter image description here

In other cases the area of the selected boundary was the outer:

enter image description here

I do know that there is the left-hand convention for storing inner area of a polygon in sql server as spatial data. Clearly fixing each polygon one by one will be a waste of time.

GeoJson files aren't bounded to the left-hand rule, so I need a way to fix the coordinates orientation of the polygon before storing them in the db.

I want the orientation fix logic to be done in the c# code (MVC) not by the sql server.

I've read about a possible solution for convex polygons that handels the issue with cross-product, but still this doesn't work.

  1. How can I solve this issue?
  2. Is there a program that fixes the geoJson orientation of the polygons?
  3. Is there an algorithm based on a mathematical principle?
  4. Without going to the fuss of checking the orientation and all the extra coding, is there a site where I can download world data coordinates with all the administrative levels that use the left-hand convention for polygons?

Code used to get cross product

public decimal GetCrossProductZDirectionAndLength(VectorCoordinate coordA, VectorCoordinate coordB, VectorCoordinate coordC) 
{
    //***********************************************
    //This method gets the cross product 
    //for vectors:
    //BA X BC
    //***********************************************
    decimal crossProductDirectionLength = 0;

    //These vector coordinates are relative to the origin (0,0,0)
    VectorCoordinate vectorCoordinateBA = new VectorCoordinate();
    VectorCoordinate vectorCoordinateBC = new VectorCoordinate();

    vectorCoordinateBA.X = coordA.X - coordB.X;
    vectorCoordinateBA.Y = coordA.Y - coordB.Y;
    vectorCoordinateBA.Z = coordA.Z - coordB.Z;

    vectorCoordinateBC.X = coordC.X - coordB.X;
    vectorCoordinateBC.Y = coordC.Y - coordB.Y;
    vectorCoordinateBC.Z = coordC.Z - coordB.Z;

    crossProductDirectionLength = GetDeterminat2x2(vectorCoordinateBA.Y, vectorCoordinateBA.Z, vectorCoordinateBC.Y, vectorCoordinateBC.Z);
    crossProductDirectionLength += GetDeterminat2x2(vectorCoordinateBA.X, vectorCoordinateBA.Z, vectorCoordinateBC.X, vectorCoordinateBC.Z);
    crossProductDirectionLength += GetDeterminat2x2(vectorCoordinateBA.X, vectorCoordinateBA.Y, vectorCoordinateBC.X, vectorCoordinateBC.Y);

    return crossProductDirectionLength;
}


private decimal GetDeterminat2x2(decimal a, decimal b, decimal c, decimal d) 
{
    //get the determinat for a 2x2 matrix
    //|a b|
    //|c d|

    return (a*d)-(b*c);
}
  • Could you expand on the solution you have already tried and what about it did not work? – MaryBeth Nov 5 '15 at 18:00
  • @MaryBeth I've included the code and an answer, please check if it makes sense, I've tested it on a few polygons and it seems to work. – Luther Nov 6 '15 at 2:01
2

I don't have a C# solution but here is how to do it in Java using JTS and GeoTools. But you should be able to recreate it in any language which provides some basic libraries/methods.

The algorithm comes down to

for each polygon do
   if outer ring is counter clockwise then
      reverse outer ring 
      for each inner ring
         reverse it

so in java

while (it.hasNext()) {
    SimpleFeature f = (SimpleFeature) it.next();
    Geometry geom = (Geometry) f.getDefaultGeometry();
    System.out.println(geom);
    if (geom instanceof Polygon) {
        f.setDefaultGeometry(fixPolygon((Polygon) geom));
    } else if (geom instanceof MultiPolygon) {
        MultiPolygon multi = (MultiPolygon) geom;
        int numGeometries = multi.getNumGeometries();
        Polygon[] polys = new Polygon[numGeometries];
        for (int i = 0; i < numGeometries; i++) {
            polys[i] = fixPolygon((Polygon) multi.getGeometryN(i));
        }
        f.setDefaultGeometry(GEOMFAC.createMultiPolygon(polys));
    }
    ret.add(f);
}


private Polygon fixPolygon(Polygon geom, boolean cw) {
    LineString ring = geom.getExteriorRing();
    LinearRing extRing;
    if (RobustCGAlgorithms.isCCW(ring.getCoordinates()) == cw) {
        extRing = JTSUtilities.reverseRing((LinearRing) ring);
    } else {
        extRing = (LinearRing) ring;
    }
    Polygon ret;
    int numInteriorRing = geom.getNumInteriorRing();
    if (numInteriorRing > 0) {
        LinearRing[] holes = new LinearRing[numInteriorRing];
        for (int i = 0; i < numInteriorRing; i++) {
            LineString inner = geom.getInteriorRingN(i);
            if (RobustCGAlgorithms.isCCW(inner.getCoordinates()) != cw) {
                holes[i] = JTSUtilities.reverseRing((LinearRing) inner);
            } else {
                holes[i] = (LinearRing) inner;
            }
        }
        ret = GEOMFAC.createPolygon(extRing, holes);
    } else {
        ret = GEOMFAC.createPolygon(extRing);
    }
    return ret;
}
1

I've tried a possible solution to the convex/concave polygon coordinate orientation, based on the left-hand rule (right-hand if required).

Basically we can find the area of a polygon by adding the areas of the trapezoids defined by the polygons edge and a line corresponding to the min Y value of all the coordinates of the polygon. The min Y line is parallel to the X axis.

enter image description here

The area of the trapezoid can be found by applying and rearranging the area formula:

area = [(x2 - x1) * (y2 + y1)] / 2

enter image description here

When the program adds up all of the trapezoid areas, the sides on the polygon’s bottom give negative areas because x1 > x2 (this depends on the coordinates orientation). Those areas cancel out the parts of the other trapezoids that lie outside of the polygon. This method gives strange results for self-intersecting polygons.

The code then loops over the polygon’s segments, calculates the area under each, adds them up, and returns the total.

The total calculated area is negative if the polygon is oriented clockwise.

Clockwise or counter-clockwise check can be easily achieved with this solution.


  public class VectorCoordinate
{
    public decimal X { get; set; }
    public decimal Y { get; set; }
}


public class VectorHelper 
{
    public decimal GetAreaOfPolygon(List<VectorCoordinate> listPolygonCoordinate) 
    {
        decimal minYcoord = GetMinYCoordFromPolygon(listPolygonCoordinate);
        //Get the first coord to check if the last coord is the same
        //Area of the polygon
        decimal areaOfPolygon = 0;
        //List of corrds to array
        VectorCoordinate[] vectorCoordinate = listPolygonCoordinate.ToArray(); 

        for (int i = 0; i < vectorCoordinate.Length-1; i++)
        {
            areaOfPolygon += GetAreaOfTrapezoid(vectorCoordinate[i].X, vectorCoordinate[i].Y, vectorCoordinate[i + 1].X, vectorCoordinate[i + 1].Y, minYcoord);
        }
        //posite/negative area result will tell us if the orientation is clockwise or counterclockwise 
        return areaOfPolygon;
    }

    private decimal GetMinYCoordFromPolygon(List<VectorCoordinate> listPolygonCoordinate) 
    {
        decimal minYcoord = listPolygonCoordinate.First().Y;

        foreach(var polygonCoord in listPolygonCoordinate)
        {
            if (polygonCoord.Y < minYcoord)
                minYcoord = polygonCoord.Y;
        }

        return minYcoord;
    }
    private decimal GetAreaOfTrapezoid(decimal x1, decimal y1, decimal x2, decimal y2, decimal minYcoord) 
    {
        return ((x2 - x1) * ( (y2-minYcoord) + (y1-minYcoord))) / 2;
    }
}

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.