Pages: 1
- Sujet précédent - Recherche moyen de sélectionner polygones selon leurs orientations - Sujet suivant
#1 Wed 14 April 2010 10:11
- b.ferry
- Juste Inscrit !
- Date d'inscription: 14 Apr 2010
- Messages: 5
Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour,
dans le cadre d'une étude sur les potentialités solaires en toiture au sein d'une communauté d'agglomération, nous recherchons le moyen de sélectionner seulement certains polygones (qui représentent les batiments) de la base topo.
Le but est de ne conserver que les batiments orientés dans l'axe Nord-Sud.
Il semble qu'une simple requête SQL ne soit pas suffisante pour réaliser cela. Quelqu'un voit-il un moyen de faire ce genre de requête ? Par exemple en admettant que l'on prend la longueur max d'un polygone rectangulaire d'orientation Est-Ouest cela donne une toiture inclinée face Nord et Sud (ce qui nous interesse) mais comment faire ce genre de requête ?
Je vous remercie d'avance pour votre aide.
Cordialement.
Benoît.
Hors ligne
#2 Wed 14 April 2010 10:58
Re: Recherche moyen de sélectionner polygones selon leurs orientations
bonjour,
quel logiciel utilisez-vous?
cordialement
Xavier Germain
SARL Géodésie: Cartographie Numérique
2 Résidence de Hent Du
29170 Fouesnant
Hors ligne
#3 Wed 14 April 2010 11:15
- b.ferry
- Juste Inscrit !
- Date d'inscription: 14 Apr 2010
- Messages: 5
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Le logiciel Mapinfo V.9.x
Hors ligne
#4 Wed 14 April 2010 13:35
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour,
Un idée qui vaut ce qu'elle vaut, sans rapport avec MapInfo.
Si l'on prend le rectangle encadrant du polygone(xmin,ymin,xmax,ymax) le vecteur [(xmin,ymin);(xmax,ymax)] va vous donner l'orientation générale selon la plus grande longueur, non ?
A+
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#5 Wed 14 April 2010 15:22
- Jean-Jacques
- Participant actif
- Lieu: Aurillac
- Date d'inscription: 4 Jan 2006
- Messages: 99
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour,
Avec mes moyens:
Faire des courbes iso-orientation : MapInfo V10 / BD Alti / Vertical Mapper
puis Requête SQL entre Bâti et Courbes transformées en polygones selon l'orientation désirée : Bati de la BDTopo ou Parcellaire
encore faut-il avoir des infos sur les pans de toitures (un toit à 45 ° tourné au nord ne verra pas le soleil longtemps même si
le bâtiment est plein sud)
Cordialement
Jean-Jacques
DDT15
Jean-Jacques
DDT15
Hors ligne
#6 Wed 14 April 2010 19:14
- Damien BEAUSEIGNEUR
- Participant assidu
- Lieu: meyzieu
- Date d'inscription: 5 Sep 2005
- Messages: 425
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour
christophe désolé mais ça ne marche pas, ou ça va fausser le résultat.
l'angle donné pour un rectangle orienté pour sa longueur dans le sens Est-Ouest donnera un angle supérieur à 0°.
et pour un carré et ceci quelque soit sa position donnera un angle de 45°.
La question est quelle précision angulaire voulez-vous?
l'idée de l'orientation sera donnée par la différence xmax-xmin par rapport à
si xmax-xmin > ymax-ymin alors l'angle est de moins de 45°
si xmax-xmin < ymax-ymin alors l'angle est de plus de 45°
pour être plus fin il faut travailler avec tous les segments du contour du polygone
en calculant la longueur la plus grande entre les segments étant orienté avec un angle supérieur à 45°
et ceux orienté avec un angle inférieur à 45°
voila pour l'idée.
ensuite effectivement avec certains pans de toits on peut avoir 4 faces ou voir pire pratiquement uniquement orienté plein nord, ce qui n'est pas très pratique pour des panneaux solaires.
Pour tout ça. l'aide d'une photo aérienne peut aider.
cordialement.
Hors ligne
#7 Thu 15 April 2010 08:16
- ChristopheV
- Membre
- Lieu: Ajaccio
- Date d'inscription: 7 Sep 2005
- Messages: 3175
- Site web
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour,
christophe désolé mais ça ne marche pas, ou ça va fausser le résultat.
Il n'y a pas à être désolé Effectivement ça va pas marcher.
A+
Christophe
L'avantage d'être une île c'est d'être une terre topologiquement close
Hors ligne
#8 Fri 16 April 2010 15:59
- b.ferry
- Juste Inscrit !
- Date d'inscription: 14 Apr 2010
- Messages: 5
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Bonjour
christophe désolé mais ça ne marche pas, ou ça va fausser le résultat.
l'angle donné pour un rectangle orienté pour sa longueur dans le sens Est-Ouest donnera un angle supérieur à 0°.
et pour un carré et ceci quelque soit sa position donnera un angle de 45°.
La question est quelle précision angulaire voulez-vous?
l'idée de l'orientation sera donnée par la différence xmax-xmin par rapport à
si xmax-xmin > ymax-ymin alors l'angle est de moins de 45°
si xmax-xmin < ymax-ymin alors l'angle est de plus de 45°
pour être plus fin il faut travailler avec tous les segments du contour du polygone
en calculant la longueur la plus grande entre les segments étant orienté avec un angle supérieur à 45°
et ceux orienté avec un angle inférieur à 45°
voila pour l'idée.
ensuite effectivement avec certains pans de toits on peut avoir 4 faces ou voir pire pratiquement uniquement orienté plein nord, ce qui n'est pas très pratique pour des panneaux solaires.
Pour tout ça. l'aide d'une photo aérienne peut aider.
cordialement.
Votre explication est très intéressante et me permet d'y voir plus clair donc je vous en remercie.
Cependant j'ai un niveau débutant en SIG, donc j'ai encore plein d'autres questions
Donc visiblement xmin xmax ymin ymax correspondent aux coordonnées de mon polygone. J'ai fait des recherches et visiblement il existe des moyens de les extraires :
- la fonction Objectgeography() permet de récupérer les coordonnées d'un polygone.
- la fonction MBR() permet de récupérer un object rectangle incluant mon polygone.
En utilisant ces deux fonctions, il y a moyen de récupérer les coordonnées de mes polygones et prendre ceux qui sont compris entre 0 et 45° (0 = horizontal ; 45 ° = orientation maximale intéréssante)
Le soucis c'est que je n'ai jamais utilisé mapbasic. Je ne sais donc pas programmer un script qui me permettrait de réaliser çà. Si jamais ce n'est pas très compliqué et trop long , quelqu'un sait-il quel code je dois taper ? Les 2 personnes du service SIG de ma collectivité n'ont pas le temps de m'aider malheureusement...
Merci à tous ceux qui ont déjà repondu.
Benoît
Hors ligne
#9 Tue 20 April 2010 09:44
- b.ferry
- Juste Inscrit !
- Date d'inscription: 14 Apr 2010
- Messages: 5
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Ce sujet n'inspire plus personne ?
J'ai trouvé un code déjà tout fait pour réaliser ce genre de requête sous Arcview malheureusement pas sous Mapbasic
Hors ligne
#10 Tue 20 April 2010 11:48
- Spacejo
- Membre
- Lieu: Nancy
- Date d'inscription: 17 Aug 2008
- Messages: 2511
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Salut,
serait 'il possible d'avoir un lien vers ce code ou peux tu nous le poster?
On peut peut être l'adapter pour MI
A+
Joël
Hors ligne
#11 Tue 20 April 2010 16:03
- b.ferry
- Juste Inscrit !
- Date d'inscription: 14 Apr 2010
- Messages: 5
Re: Recherche moyen de sélectionner polygones selon leurs orientations
Ok, tout d'abord voici le lien du script : http://arcscripts.esri.com/details.asp?dbid=14570
Summary
The purpose of the Polygon Diameter Azimuth Tool is to find the azimuth of a polygons diameter (major-axis, greatest antipodal distance…) and write the azimuth value to an ‘Azimuth’ field. An option has been added that creates a ‘Cardinal’ field and populates this field with a 90 Degree rotation based upon whether the Height of a polygon’s envelope is greater than its width. If a polygons envelope is wider than it is high then the ‘Cardinal’ field is set to a rotation of 90 degrees. All of the above is done to be able to rotate map pages using ESRI’s sample Map Book application. The Polygon Diameter Azimuth Tool requires ArcMap ArcView. This tool runs on a polygon layer. The Tool exists as a VBA form that needs to be imported in the VBA Editor. The form is called “frm_FindDiameter”. This tool was already well developed (better than I have done) by William Huber and Jeff Jenness in ArcView 3.X. I am just standing on their shoulders.
Il s'agit exactement de ce que je souhaiterais pouvoir faire.
Le code est un peu long mais le voici :
Code:
VERSION 5.00 Begin {C62A69F0-16DC-11CE-9E98-00AA00574A4F} frm_FindDiameter Caption = "Find Polgon Diameters" ClientHeight = 3120 ClientLeft = 45 ClientTop = 330 ClientWidth = 4710 OleObjectBlob = "frm_FindDiameter.frx":0000 StartUpPosition = 1 'CenterOwner End Attribute VB_Name = "frm_FindDiameter" Attribute VB_GlobalNameSpace = False Attribute VB_Creatable = False Attribute VB_PredeclaredId = True Attribute VB_Exposed = False Public turncomplete As Boolean Private Sub cbx_Sel_Layer_Change() Dim pMxDoc As IMxDocument Dim pLayer As ILayer Dim pFeatLayer As IFeatureLayer Dim pEnumLayers As IEnumLayer Dim pFeatClass As IFeatureClass Dim strLayer As String Dim pAV As IActiveView Dim pFields As IFields Dim pField As IField Dim pEnumField As IEnumFieldMap Dim myArr As IArray chkbx_CreateAz.Enabled = True chkbx_CreateCard.Enabled = True Set pMxDoc = Application.Document Set pAV = pMxDoc.FocusMap Set pEnumLayers = pMxDoc.FocusMap.Layers pEnumLayers.Reset Set pLayer = pEnumLayers.Next Do While Not pLayer Is Nothing If UCase(pLayer.Name) = UCase(cbx_Sel_Layer.Value) Then Exit Do Set pLayer = pEnumLayers.Next Loop If pLayer Is Nothing Then Exit Sub End If Set pFeatLayer = pLayer Set pFeatClass = pFeatLayer.FeatureClass IndexA = pFeatClass.FindField("Azimuth") If (IndexA < 0) Then chkbx_CreateAz.Enabled = True End If IndexA = pFeatClass.FindField("Cardinal") If (IndexA < 0) Then chkbx_CreateCard.Enabled = True End If End Sub Private Sub chkbx_CreateAz_Click() If chkbx_CreateAz Then btn_OK.Enabled = True Else If chkbx_CreateCard = False Then btn_OK.Enabled = False End If End If End Sub Private Sub chkbx_CreateCard_Click() If chkbx_CreateCard Then btn_OK.Enabled = True Else If chkbx_CreateAz = False Then btn_OK.Enabled = False End If End If End Sub Private Sub UserForm_Initialize() Dim pDoc As IMxDocument Dim pELayer As IEnumLayer Dim pLayer As ILayer, pLayer2 As ILayer Dim pFLayer As IFeatureLayer Dim pFClass As IFeatureClass Dim pMap As IMap Dim i As Integer Set pDoc = ThisDocument Set pELayer = pDoc.ActiveView.FocusMap.Layers pELayer.Reset Set pLayer = pELayer.Next Do While Not pLayer Is Nothing If TypeOf pLayer Is IFeatureLayer Then Set pFLayer = pLayer If pFLayer.FeatureClass.ShapeType = esriGeometryPolygon Then cbx_Sel_Layer.AddItem pLayer.Name End If End If Set pLayer = pELayer.Next Loop If cbx_Sel_Layer.ListCount = 0 Then MsgBox "Tool needs polygon layer", vbCritical, "EXITING" Exit Sub End If ' Set controls on Form chkbx_CreateAz.Enabled = False chkbx_CreateCard.Enabled = False btn_OK.Enabled = False End Sub Private Sub btn_OK_Click() Dim strLayer As String Dim strFile As String strLayer = cbx_Sel_Layer.SelText ' Set controls on Form Call ProcessSelections(strLayer) MsgBox "Done", vbOKOnly, "Successfully Completed" 'frm_FindDiameter.Hide End Sub Public Sub ProcessSelections(sLayerName As String) Dim pMxDoc As IMxDocument Dim pLayer As ILayer Dim pFeatLayer As IFeatureLayer Dim pEnumLayers As IEnumLayer Dim pFeatClass As IFeatureClass Dim pFeatCursor As IFeatureCursor Dim pFeat As IFeature Dim nFeat As Integer Dim pTopoOp As ITopologicalOperator Dim pHull As IPolygon Dim pPolygon As IPolygon Dim pPtColl As IPointCollection Dim PolyDiPts As IPointCollection Dim PolyDiMP As IMultipoint Dim pAV As IActiveView Dim pGC As IGraphicsContainer Dim pline As ILine Dim getRadian As Double Dim getAngle As Integer, iAriAngle As Integer, iGeoAngle As Integer Dim pi As Double Dim IndexA As Long Dim pField As IFieldEdit pi = 4 * Atn(1) 'Get Calculated Value of pi Set pMxDoc = Application.Document Set pAV = pMxDoc.FocusMap Set pGC = pAV Set pEnumLayers = pMxDoc.FocusMap.Layers Set PolyDiMP = New Multipoint Set PolyDiPts = PolyDiMP ' find the requested layer: Set pLayer = pEnumLayers.Next Do While Not pLayer Is Nothing If UCase(pLayer.Name) = UCase(sLayerName) Then Exit Do Set pLayer = pEnumLayers.Next Loop If pLayer Is Nothing Then Exit Sub End If ' get the feature cursor: Set pFeatLayer = pLayer Set pFeatClass = pFeatLayer.FeatureClass If chkbx_CreateAz Then IndexA = pFeatClass.FindField("Azimuth") If IndexA < 0 Then Set pField = New Field With pField .Type = esriFieldTypeInteger .Name = "Azimuth" End With pFeatClass.AddField pField End If Index_Az = pFeatClass.FindField("Azimuth") End If If chkbx_CreateCard Then IndexA = pFeatClass.FindField("Cardinal") If IndexA < 0 Then Set pField = New Field With pField .Type = esriFieldTypeInteger .Name = "Cardinal" End With pFeatClass.AddField pField End If Index_Card = pFeatClass.FindField("Cardinal") End If Set pFeatCursor = pFeatClass.Search(Nothing, False) 'Get a Cursor for all polygons in layer ' loop through the features: Set pFeat = pFeatCursor.NextFeature 'Get first polygon If (pFeat.Shape.GeometryType = esriGeometryPolygon) Then 'make sue it is a polgyon nFeat = 0 'Initialize polygon counter Do While Not pFeat Is Nothing 'Loop through all polygons nFeat = nFeat + 1 'Keep track of number of polygon processed Set pFeatPoints = pFeat.Shape 'Get the points that make up a polygon Set pPolygon = pFeat.ShapeCopy Set pTopoOp = pFeatPoints 'Assign points to a Topological Operator which allows for constructing new geometries Set pHull = pTopoOp.ConvexHull 'Get the convex hull of the point collection and assign to a polygon Set PolyDiPts = GetPolyDiPts(pHull) 'Pass the convex hull polgyon to a function to retrieve antipodal points of maximum distance apart Dim Apt As IPoint Dim Bpt As IPoint Set Apt = New Point Set Bpt = New Point Set pline = New Line 'Create empty line object PolyDiPts.QueryPoint 0, Apt 'Assign first antipodal point to a new point PolyDiPts.QueryPoint 1, Bpt 'Assgin second antipodal point to a new Point pline.PutCoords Apt, Bpt 'Create line from antipodal points getRadian = pline.Angle 'Retrieve angle of line object from a property of that object If getRadian < 0 Then 'Convert radians to Polar iAriAngle = 360 - (((getRadian * 360) * -1) / (2 * pi)) Else iAriAngle = (getRadian * 360) / (2 * pi) End If getAngle = 45 - (iAriAngle - 45) If getAngle < 0 Then 'Convert polar to Azimuth iGeoAngle = getAngle + 360 Else iGeoAngle = getAngle End If If chkbx_CreateAz Then pFeat.Value(Index_Az) = iGeoAngle 'Set Attribute field of polygon to Azimuth pFeat.Store End If If chkbx_CreateCard Then pFeat.Value(Index_Card) = GetCard(pPolygon) 'Set Attribute field of polygon to Azimuth pFeat.Store End If Set pFeat = pFeatCursor.NextFeature 'Select next polygon to process Loop Else MsgBox ("Wrong type feature geometry") Exit Sub End If End Sub Private Function GetCard(pPolygon As IPolygon) As Integer 'Get envelope Dim pArea As IArea Dim myCntrd As IPoint Dim ratio As Double Dim pEnv As IEnvelope Set pArea = pPolygon Set myCntrd = pArea.Centroid Set pEnv = pPolygon.Envelope ratio = pEnv.Height / pEnv.Width If ratio > 1 Then GetCard = 0 Else GetCard = 270 End If End Function Private Function MakeMarkerElement(pPoint As IPoint) As IMarkerElement Dim pElement As IElement Dim pMxDoc As IMxDocument Dim pDocDefaultSymbols As IDocumentDefaultSymbols Set pElement = New MarkerElement pElement.Geometry = pPoint Set MakeMarkerElement = pElement Set pMxDoc = ThisDocument Set pDocDefaultSymbols = pMxDoc MakeMarkerElement.Symbol = pDocDefaultSymbols.MarkerSymbol Exit Function End Function Public Function GetPolyDiPts(pHull As IPolygon) As IPointCollection 'Passed convex hull polygon. returns max antipodal pair Dim pPtColl As IPointCollection Dim AntiPodalNdex As Integer Dim maxDist As Double Dim myPtColl As IPointCollection Dim MaxDiColl As IPointCollection Dim MaxMP As IMultipoint Set MaxDiColl = MaxMP 'Set to empty collection of multipoint object Set myPtColl = MaxMP 'Set to empty collection of multipoint object Set GetPolyDiPts = MaxMP 'Set to empty collection of multipoint object Set pPtColl = pHull 'Get points that make up convex hull polygon pPtColl.RemovePoints pPtColl.PointCount - 1, 1 'Remove repeated polygon closure point If (pPtColl.PointCount = 3) Then 'If poygon is a triangle Set GetPolyDiPts = FindPtsFor3(pPtColl, 0, 1, 2) 'Call function that solves easy case. Pass points and index starting points. Returns antipodal pair. Else 'If polygon is not triangle Dim over As Boolean 'Need to define Dim currentIndex As Integer 'Algorithm needs to track 3 points in point collection Dim prevIndex As Integer 'previous index number to current point in point collection Dim Dist As Long 'Need to define over = False currentIndex = 0 'Initialize current point index in point collection prevIndex = FindPrvCCWpt(pPtColl, currentIndex) 'Find the previous index value of current index point. Passes point collection and current index point. returns previous point index value AntiPodalNdex = FindAntiPndx(pPtColl, currentIndex, 1) 'Find the index value to the antipodal point match for the current index point maxDist = FindMaxFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'Finds maxDist of Antipodal pair (pair could be previous index or current index and Antpodal index Set MaxDiColl = FindPtsFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'Determine antipodal pair from two posibilities Dist = 0 'Initialize max distance While (over = False) 'Loop through other points in collection traversing convex hull prevIndex = currentIndex 'increment previous index currentIndex = FindNxtCCWpt(pPtColl, currentIndex) 'increment current index Dist = FindMaxFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'We can now start at previous antipodal index value when determining next antipodal point Set myPtColl = FindPtsFor3(pPtColl, currentIndex, prevIndex, AntiPodalNdex) 'Get antipodal pair If (Dist > maxDist) Then 'Compare distances of antipodal pairs maxDist = Dist 'Assign max distance if greater Set MaxDiColl = myPtColl 'Assign point collection with max Antipodal pair End If If turncomplete Then over = True 'Assign boolean to indicate traversing convex hull is complete Wend Set GetPolyDiPts = MaxDiColl 'return point collection of antipodal points representing polygon diameter End If End Function Public Function Euclidean_Dist(Firstpt As IPoint, Secondpt As IPoint) As Double 'Passed points calculates distance. Returns distance Euclidean_Dist = (Sqr(((Firstpt.X - Secondpt.X) ^ 2) + (Firstpt.Y - Secondpt.Y) ^ 2)) 'Pythagorean Theorem End Function Public Function FindPtsFor3(pPtColl As IPointCollection, ndex1 As Integer, ndex2 As Integer, ndex3 As Integer) As IPointCollection 'Passed polygon point collection and index values of points. Returns antipodal pair Dim dist_12 As Double Dim dist_23 As Double Dim dist_31 As Double Dim pt1 As IPoint Dim pt2 As IPoint Dim pt3 As IPoint Dim Chk As Integer Dim myMP As IMultipoint Dim myPtColl As IPointCollection Set myMP = New Multipoint Set myPtColl = myMP Set pt1 = New esriGeometry.Point Set pt2 = New esriGeometry.Point Set pt3 = New esriGeometry.Point pPtColl.QueryPoint ndex1, pt1 'Get point in point collection from passed index value pPtColl.QueryPoint ndex2, pt2 'Get point in point collection from passed index value pPtColl.QueryPoint ndex3, pt3 'Get point in point collection from passed index value dist_12 = Euclidean_Dist(pt1, pt2) 'Call function to get distance between points dist_23 = Euclidean_Dist(pt2, pt3) 'Call function to get distance between points dist_31 = Euclidean_Dist(pt3, pt1) 'Call function to get distance between points Dim maxDist As Double maxDist = dist_12 'Assign initial maximum distance Chk = 1 If (dist_23 > maxDist) Then 'Check to see which distance is greatest Chk = 2 ElseIf (dist_31 > maxDist) Then Chk = 3 End If Select Case Chk 'Depending on Case assign max antipodal pair to point collection Case 1 myPtColl.AddPoint pt1 myPtColl.AddPoint pt2 Case 2 myPtColl.AddPoint pt2 myPtColl.AddPoint pt3 Case 3 myPtColl.AddPoint pt3 myPtColl.AddPoint pt1 End Select Set FindPtsFor3 = myPtColl 'Assign point collection to returning value End Function Public Function FindMaxFor3(pPtColl As IPointCollection, ndex1 As Integer, ndex2 As Integer, ndex3 As Integer) As Double ' Dim dist_23 As Double Dim dist_31 As Double Dim pt1 As IPoint Dim pt2 As IPoint Dim pt3 As IPoint Dim Chk As Integer Set pt1 = New esriGeometry.Point Set pt2 = New esriGeometry.Point Set pt3 = New esriGeometry.Point pPtColl.QueryPoint ndex1, pt1 'Get points from their index position pPtColl.QueryPoint ndex2, pt2 pPtColl.QueryPoint ndex3, pt3 dist_12 = Euclidean_Dist(pt1, pt2) 'Calculate distance between points dist_23 = Euclidean_Dist(pt2, pt3) dist_31 = Euclidean_Dist(pt3, pt1) Dim maxDist As Double maxDist = dist_12 'Initialize maxDist value If (dist_23 > maxDist) Then maxDist = dist_23 'Determine maximum distance If (dist_31 > maxDist) Then maxDist = dist_31 FindMaxFor3 = maxDist 'Return maximum distance End Function Public Function FindNxtCCWpt(pPtColl As IPointCollection, ndex As Integer) As Integer 'passes point collection and point index value ' Find Next Counter ClockWise Point If (ndex >= pPtColl.PointCount - 1) Then 'If finished traversing convex hull turncomplete = True 'Notify calling routines that traversing is complete FindNxtCCWpt = 0 Else FindNxtCCWpt = ndex + 1 End If End Function Public Function FindPrvCCWpt(pPtColl As IPointCollection, ndex As Integer) As Integer 'Passed point collection and current point index value. Returns previous point index value 'Find Previous Count Clockwise Point FindPrvCCWpt If (ndex <= 0) Then 'If at the 0 index point then FindPrvCCWpt = pPtColl.PointCount - 1 'Get last point in point collection index value Else 'if index not 0 FindPrvCCWpt = ndex - 1 'Get the previous point index value End If End Function Public Function FindAntiPndx(pPtColl As IPointCollection, ndex As Integer, startndex As Integer) As Integer 'Passed the point collection, current point index value, index value to start process to find antipodal point Dim a As Integer Dim b As Integer Dim c As Integer Dim current As Integer Dim area_abc As Double Dim area_abcurrent As Double Dim Apt As IPoint 'Establish point objects to work with Dim Bpt As IPoint Dim Cpt As IPoint Set Apt = New esriGeometry.Point Set Bpt = New esriGeometry.Point Set Cpt = New esriGeometry.Point b = ndex 'Value of current point index c = startndex 'Value of starting index point for processing a = FindPrvCCWpt(pPtColl, b) 'Find the previous point index value of current index current = FindPrvCCWpt(pPtColl, c) 'find the previous point index of the starting index for processing pPtColl.QueryPoint a, Apt 'Get the points from the point collection for processing pPtColl.QueryPoint b, Bpt pPtColl.QueryPoint c, Cpt area_abc = AreaOfTriangle(Apt, Bpt, Cpt) 'Calculate the area the intial 3 points pPtColl.QueryPoint current, Cpt 'Get index starting point area_abcurrent = AreaOfTriangle(Apt, Bpt, Cpt) 'Calculate area While (area_abc >= area_abcurrent) 'While traversing the points of a convex hull polygon current = c 'calculating area from current and previous point and next traversed point c = FindNxtCCWpt(pPtColl, c) 'When the area of the triangle starts decreasing you have found your pPtColl.QueryPoint c, Cpt 'antipodal pair of points area_abc = AreaOfTriangle(Apt, Bpt, Cpt) pPtColl.QueryPoint current, Cpt area_abcurrent = AreaOfTriangle(Apt, Bpt, Cpt) Wend FindAntiPndx = current 'Return the index value of the antipodal point of the point collection End Function Public Function AreaOfTriangle(Apt As IPoint, Bpt As IPoint, Cpt As IPoint) As Double 'Pass three points return area 'Calculate area of triangle using the Heron's formula Dim a As Double Dim b As Double Dim c As Double Dim s As Double a = Abs(Euclidean_Dist(Apt, Bpt)) 'Get distance between points b = Abs(Euclidean_Dist(Bpt, Cpt)) c = Abs(Euclidean_Dist(Cpt, Apt)) s = 0.5 * (a + b + c) ' Semi-perimeter AreaOfTriangle = Sqr(s * ((s - a) * (s - b) * (s - c))) End Function Public Function Features2ConvexHull(sLayerName As String, sTextFileName As String, Optional nMaxFeatures As Double) Dim pFeat As IFeature Dim pLayer As ILayer Dim pMxDoc As IMxDocument Dim pEnumLayers As IEnumLayer Dim i As Integer Dim pFeatClass As IFeatureClass Dim pFeatLayer As IFeatureLayer Dim pFeatCursor As IFeatureCursor Dim sY As String, sX As String, sFeatN As String, sFeatVerts As String Dim strCenX As String, strCenY As String, strAngle As String Dim lFileID As Long Dim pFeatPoints As IPointCollection Dim pArea As IArea Dim pPoints As Point Dim nFeat As Double Dim pAV As IActiveView Dim pGC As IGraphicsContainer ' get the document and map objects Set pMxDoc = Application.Document Set pAV = pMxDoc.FocusMap Set pGC = pAV Set pEnumLayers = pMxDoc.FocusMap.Layers ' find the requested layer: Set pLayer = pEnumLayers.Next Do While Not pLayer Is Nothing If UCase(pLayer.Name) = UCase(sLayerName) Then Exit Do Set pLayer = pEnumLayers.Next Loop If pLayer Is Nothing Then Exit Function End If ' get the feature cursor: Set pFeatLayer = pLayer Set pFeatClass = pFeatLayer.FeatureClass Set pFeatCursor = pFeatClass.Search(Nothing, False) ' loop through the features: Set pFeat = pFeatCursor.NextFeature If (pFeat.Shape.GeometryType = esriGeometryPolygon) Then Do While Not pFeat Is Nothing nFeat = nFeat + 1 If (nFeat > nMaxFeatures) And (nMaxFeatures > 0) Then Exit Do End If Set pFeatPoints = pFeat.Shape Dim pTopoOp As ITopologicalOperator Set pTopoOp = pFeatPoints Dim pHull As IPolygon Set pHull = pTopoOp.ConvexHull Loop Else MsgBox ("Wrong type feature geometry") Return End If End Function
Ca me parait assez long pour pouvoir être adapté sous mapinfo, peut etre (et même surement) qu'il existe déjà le même script sous mapbasic mais je ne l'ai pas trouvé.
Merci pour votre aide en tout cas.
Hors ligne
Pages: 1
- Sujet précédent - Recherche moyen de sélectionner polygones selon leurs orientations - Sujet suivant