Suite

GeoJson to Spatial Geography dans SQL Server 2012, correction de l'orientation du polygone

GeoJson to Spatial Geography dans SQL Server 2012, correction de l'orientation du polygone


J'ai téléchargé depuis Mapzen de nombreux fichiers geoJson liés aux pays avec leur niveau d'administration respectif. Ensuite, j'ai analysé tous les fichiers geoJson et les ai enregistrés sur mon serveur SQL 2012 en tant que données spatiales géographiques.

J'ai ensuite vérifié la carte du Luxembourg et j'ai remarqué des résultats étranges. Certaines des cartes ont été créées correctement, avec la zone intérieure correspondante :

Dans d'autres cas, l'aire de la limite sélectionnée était l'extérieur :

Je sais qu'il existe la convention de gauche pour stocker la zone intérieure d'un polygone dans le serveur SQL en tant que données spatiales. Il est clair que fixer chaque polygone un par un sera une perte de temps.

Les fichiers GeoJson ne sont pas limités à la règle de gauche, j'ai donc besoin d'un moyen de fixer l'orientation des coordonnées du polygone avant de les stocker dans la base de données.

Je souhaite que la logique de correction d'orientation soit effectuée dans le code c# (MVC) et non par le serveur SQL.

J'ai lu une solution possible pour les polygones convexes qui résout le problème avec le produit croisé, mais cela ne fonctionne toujours pas.

  1. Comment puis-je résoudre ce problème ?
  2. Existe-t-il un programme qui fixe l'orientation geoJson des polygones ?
  3. Existe-t-il un algorithme basé sur un principe mathématique ?
  4. Sans aller jusqu'à vérifier l'orientation et tout le codage supplémentaire, existe-t-il un site où je peux télécharger les coordonnées des données mondiales avec tous les niveaux administratifs qui utilisent la convention de gauche pour les polygones ?

Code utilisé pour obtenir le produit croisé

public decimal GetCrossProductZDirectionAndLength(VectorCoordinate coordA, VectorCoordinate coordB, VectorCoordinate coordC) { //************************************ ************* //Cette méthode obtient le produit croisé //pour les vecteurs : //BA X BC //***************** ****************************** décimal crossProductDirectionLength = 0; //Ces coordonnées vectorielles sont relatives à l'origine (0,0,0) VectorCoordinate vectorCoordinateBA = new VectorCoordinate(); VectorCoordinate vectorCoordinateBC = new VectorCoordinate(); vecteurCoordonnéeBA.X = coordA.X - coordB.X; vecteurCoordonnéeBA.Y = coordA.Y - coordB.Y; vecteurCoordonnéeBA.Z = coordA.Z - coordB.Z; vecteurCoordonnéeBC.X = coordC.X - coordB.X; vecteurCoordonnéeBC.Y = coordC.Y - coordB.Y; vecteurCoordonnéeBC.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(décimal a, décimal b, décimal c, décimal d) { // obtient le determinat pour une matrice 2x2 //|a b| //|c d| retourner (a*d)-(b*c); }

Je n'ai pas de solution C# mais voici comment le faire en Java en utilisant JTS et GeoTools. Mais vous devriez pouvoir le recréer dans n'importe quel langage fournissant des bibliothèques/méthodes de base.

L'algorithme se résume à

pour chaque polygone faire si l'anneau extérieur est dans le sens inverse des aiguilles d'une montre alors inverser l'anneau extérieur pour chaque anneau intérieur l'inverser

donc en java

while (it.hasNext()) { SimpleFeature f = (SimpleFeature) it.next(); Géométrie geom = (Géométrie) 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 < numGeometry; 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) anneau ; } Polygone ret ; int numInteriorRing = geom.getNumInteriorRing(); if (numInteriorRing > 0) { LinearRing[] trous = new LinearRing[numInteriorRing] ; for (int i = 0; i < numInteriorRing; i++) { LineString inner = geom.getInteriorRingN(i); if (RobustCGAlgorithms.isCCW(inner.getCoordinates()) != cw) { trous[i] = JTSUtilities.reverseRing((LinearRing) inner); } else { trous[i] = (LinearRing) inner; } } ret = GEOMFAC.createPolygon(extRing, trous); } else { ret = GEOMFAC.createPolygon(extRing); } retourne ret; }

J'ai essayé une solution possible à l'orientation des coordonnées du polygone convexe/concave, basée sur la règle de la main gauche (droite si nécessaire).

Fondamentalement, nous pouvons trouver l'aire d'un polygone en addition des aires des trapèzes défini par le bord des polygones et une ligne correspondant à la valeur Y min de toutes les coordonnées du polygone. La ligne min Y est parallèle à l'axe X.

L'aire du trapèze peut être trouvée en appliquant et en réorganisant la formule de l'aire :

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

Lorsque le programme additionne toutes les zones trapézoïdales, les côtés au bas du polygone donnent des zones négatives car x1 > x2 (cela dépend de l'orientation des coordonnées). Ces zones annulent les parties des autres trapèzes qui se trouvent à l'extérieur du polygone. Cette méthode donne des résultats étranges pour les polygones auto-sécants.

Le code parcourt ensuite les segments du polygone, calcule l'aire sous chacun, les additionne et renvoie le total.

Le total calculé surface est négatif si le polygone est orienté dans le sens des aiguilles d'une montre.

Le contrôle dans le sens des aiguilles d'une montre ou dans le sens inverse des aiguilles d'une montre peut être facilement réalisé avec cette solution.


public class VectorCoordinate { public decimal X { get; ensemble; } public décimal Y { get; ensemble; } } classe publique VectorHelper { public décimal GetAreaOfPolygon(List listPolygonCoordinate) { décimal minYcoord = GetMinYCoordFromPolygon(listPolygonCoordinate); //Obtenir le premier coord pour vérifier si le dernier coord est le même //Area of ​​the polygon decimal areaOfPolygon = 0; //Liste des corrds au tableau 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); } //Le résultat de la zone positive/négative nous dira si l'orientation est dans le sens des aiguilles d'une montre ou dans le sens inverse des aiguilles d'une montre return areaOfPolygon; } décimal privé GetMinYCoordFromPolygon(Liste listPolygonCoordinate) { décimal minYcoord = listPolygonCoordinate.First().Y; foreach(var polygonCoord in listPolygonCoordinate) { if (polygonCoord.Y < minYcoord) minYcoord = polygonCoord.Y; } return minYcoord; } privé décimal GetAreaOfTrapezoid(décimal x1, décimal y1, décimal x2, décimal y2, décimal minYcoord) { return ((x2 - x1) * ( (y2-minYcoord) + (y1-minYcoord))) / 2; } }

Voir la vidéo: ArcGIS - Spatial Join - Joins attributes based on the spatial relationship