Punkt in Polygon

[dige]

Hallo Zusammen,

da wir hier ja einige Grafikcracks haben, hoffe ich, ihr könnt mir bei folgender Aufgabenstellung helfen:

Ich habe eine Liste mit Punkten (Geokoordinaten im WGS85 Format) und mehrere Flächen (Umrisse, Polygone)
und möchte prüfen, welcher Punkt sich in welchem Polygon befindet.

Wie würdet ihr da rangehen?

Ciao Dige

[NicTheQuick]

Wir hatten sowas schon mal hier im Forum. Falls die Polygone konvex sind, wird es einfach. Falls nicht, müssen härtere Geschütze ran.

Edit: Vielleicht hilft dir ja schon das da: Polygon Include - Vereinigung,Schnitt,Differenz - Testphase

[DarkDragon]

Da es hier um geografische Daten (Länderumrisse usw.) geht vermute ich, dass es konkave, zweidimensionale Polygone sind. Wenn man es nicht allzu genau braucht, würde ich die Polygone auf ein Bild malen, jedes mit einer anderen Farbe, dann prüfen welche Farbe unter dem Punkt ist. Dann ist die Genauigkeit von der Bildauflösung abhängig.

Wenn man es genau braucht verwendet man Trapezoidal Maps, die nichts weiter sind als eine Aufteilung des Raums aller Polygone in Trapeze, indem man die Vertizen vertikal erweitert bis zur umgebenden Bounding Box.

[dige]

Danke NicTheQuick und DarkDragon für die Tipps!

Durch einen Beitrag von Infratec bin ich gerade auf die Funktion

Code: Alles auswählen

Ergebnis = IsInsidePath(x.d, y.d [, KoordinatenSystem])

aufmerksam geworden.

Wenn ich die Geokoordinaten alle in Pixelkoordinaten umrechnen würde,
wäre dann diese Funktion geeignet um so einen Test zu machen? Mit der
Vector Lib habe ich leider noch keine Erfahrung.

Vielleicht geht es auch ganz anders? Die Fläche (Polygon) hole ich mir von OpenStreetMap:
https://nominatim.openstreetmap.org/sea ... und%20GmbH

Ich will nun prüfen welche Postleitzahlen sich in dieser Fläche befinden.
Die PLZ hole ich mir hier: https://github.com/mdornseif/pyGeoDb/bl ... db_plz.txt

Ciao Dige

[alter Mann]

...ist zwar schon etwas älter...

viewtopic.php?f=8&t=18209

[ccode_new]

Hallo!

siehe:
Link von Nic

->siehe:

;+ Tests the X- and Y-coordinates whether its inside, outside or on the border of a valid polygon.
Procedure PolygonPointTest(*Polygon.Polygon, X.d, Y.d)

; Quelle: https://de.wikipedia.org/wiki/Punkt-in- ... ach_Jordan

Aber ob das in die richtige Richtung führt ?

Sind die Polygone Ländern gleichzusetzen ? Sind die Polygone auf einem eigenen Bild (z.B. Landkarte) ? Haben diese Polygone eine einheitliche Abgrenzung (auch Farbe) ?

[Michael Vogel]

Der Thread ist vielleicht nicht ganz taufrisch, aber nachdem ich auf der Suche nach einer Funktion IsPointInsidePolygon bin, bin ich auch auf diesen Beitrag gestoßen:

https://stackoverflow.com/questions/275 ... in-polygon

Solange keine größeren Regionen der Erde untersucht werden, reicht die eukldische Geometrie - denn richtig rund läuft auf der Erde ohnehin nichts


Übrigens, der folgende Code läuft auch noch nicht richtig rund...

Code: Alles auswählen

Global mousePosition.Point, mouseState, polygonPoints = 7
Global Dim polygonPoints.Point(10)


polygonPoints(0)\x = 100
polygonPoints(0)\y = 200
polygonPoints(1)\x = 350
polygonPoints(1)\y = 250
polygonPoints(2)\x = 200
polygonPoints(2)\y = 100
polygonPoints(3)\x = 550
polygonPoints(3)\y = 250
polygonPoints(4)\x = 400
polygonPoints(4)\y = 300
polygonPoints(5)\x = 350
polygonPoints(5)\y = 350
polygonPoints(6)\x = 300
polygonPoints(6)\y = 500
polygonPoints(7)\x = 250
polygonPoints(7)\y = 350


Procedure InsidePolygon(*point.point,dots,Array poly.point(1))

	; int pnpoly(int npol, float *xp, float *yp, float x, float y)
	; {
	;    int i, j, c = 0;
	;    For (i = 0, j = npol-1; i < npol; j = i++) {
	;      If ((((yp[i] <= y) && (y < yp[j])) ||
	;         ((yp[j] <= y) && (y < yp[i]))) &&
	;         (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
	;        c = !c;
	;    }
	;    Return c;
	; }

	Protected i,j,c

	While i<dots
		j=(i+dots-1)%dots
		Debug Str(i)+" "+Str(j)
			If ( ( (poly(i)\y <= *point\y) And (*point\y < poly(j)\y) ) Or
				( (poly(j)\y <= *point\y) And (*point\y < poly(i)\y) ) ) And
				( *point\x < (poly(j)\x-poly(i)\x) * (*point\y-poly(i)\y) / (poly(j)\y-poly(i)\y) + poly(i)\x)
				c!1
			EndIf
		i+1
	Wend

	ProcedureReturn c

EndProcedure

Procedure PointInPolygon(*polygon.Point, sides, *mousePosition.Point)

	Static filled
	result = #False
	

	If StartVectorDrawing(CanvasVectorOutput(0))
		
		VectorSourceColor($FFFFFFFF)
		FillVectorOutput()
		
		VectorSourceColor($FF000000)
		
		result=insidePolygon(*mousePosition,polygonPoints,polygonPoints())

		MovePathCursor(polygonPoints(0)\x,polygonPoints(0)\y)
		For i = 1 To polygonPoints-1
			AddPathLine(polygonPoints(i)\x,polygonPoints(i)\y)
		Next i
		ClosePath()
		
		If result
			FillPath()
		Else
			StrokePath(5)
		EndIf
		
		For i = 0 To polygonPoints-1
			MovePathCursor(polygonPoints(i)\x,polygonPoints(i)\y)
			DrawVectorText("("+Str(i)+") "+Str(polygonPoints(i)\x)+" | "+Str(polygonPoints(i)\y))
		Next i
		
		MovePathCursor(10,50)
		DrawVectorText(Str(*mousePosition\x)+" | "+Str(*mousePosition\y))

		FillPath()

		StopVectorDrawing()
		Delay(1)
	EndIf

	ProcedureReturn result

EndProcedure


polygonPoints(polygonPoints)=polygonPoints(0)


wFlags = #PB_Window_SystemMenu | #PB_Window_ScreenCentered
If OpenWindow(0, 0, 0, 600, 630, "PointInPolygon()", wFlags)
	CanvasGadget(0, 0, 30, 600, 600, #PB_Canvas_Container)
	TextGadget(1, 0, 0, 600, 30, "", #PB_Text_Center)

	Repeat
		event = WaitWindowEvent()
		Select event
		Case #PB_Event_CloseWindow
			appQuit = 1
		Case #PB_Event_Gadget
			If EventType() = #PB_EventType_MouseMove
				mousePosition\x = GetGadgetAttribute(0, #PB_Canvas_MouseX)
				mousePosition\y = GetGadgetAttribute(0, #PB_Canvas_MouseY)
				If mouseState <> PointInPolygon(@polygonPoints(), polygonPoints, @mousePosition)
					mouseState ! 1
					If mousestate
						txt$ = "Mouse is in..."
					Else
						txt$ = "Mouse is out..."
					EndIf
					SetGadgetText(1, Str(mousePosition\x)+" | "+Str(mousePosition\y)+" >> "+txt$)
				EndIf
			EndIf
		EndSelect
	Until appQuit
EndIf

[ccode_new]

@Michael Vogel

Übrigens, der folgende Code läuft auch noch nicht richtig rund...

Warum?

Das läuft doch schon ganz schön rund.

Code: Alles auswählen

Structure TPoint
  x.d
  y.d
EndStructure

Structure TCircle
  p.TPoint
  r.d
EndStructure


Procedure EinfacherKreis(*A.TPoint, *B.TPoint, *C.TPoint, *Kreis.TCircle)
  Protected.TPoint p1, p2
  Protected.d vB, vC, vD
  p1\x = *B\x - *A\x : p1\y = *B\y - *A\y
  p2\x = *C\x - *A\x : p2\y = *C\y - *A\y
  vB = p1\x * p1\x + p1\y * p1\y
  vC = p2\x * p2\x + p2\y * p2\y
  vD = p1\x * p2\y - p1\y * p2\x
  *Kreis\p\x = (p2\y * vB - p1\y * vC) / (2 * vD)
  *Kreis\p\y = (p1\x * vC - p2\x * vB) / (2 * vD)
  *Kreis\p\x + *A\x : *Kreis\p\y + *A\y
  *Kreis\r = Sqr(Pow(*Kreis\p\x - *A\x, 2) + Pow(*Kreis\p\y - *A\y, 2))
EndProcedure

Global mousePosition.TPoint, mouseState, polygonPoints = 360
Global Dim polygonPoints.TPoint(360)
Global.i Kreis.TCircle

polygonPoints(0)\x = 100
polygonPoints(0)\y = 200
polygonPoints(1)\x = 350
polygonPoints(1)\y = 250
polygonPoints(2)\x = 200
polygonPoints(2)\y = 100


Procedure InsidePolygon(*point.TPoint, dots, Array poly.TPoint(1))
  
  ; int pnpoly(int npol, float *xp, float *yp, float x, float y)
  ; {
  ;    int i, j, c = 0;
  ;    For (i = 0, j = npol-1; i < npol; j = i++) {
  ;      If ((((yp[i] <= y) && (y < yp[j])) ||
  ;         ((yp[j] <= y) && (y < yp[i]))) &&
  ;         (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
  ;        c = !c;
  ;    }
  ;    Return c;
  ; }
  
  Protected i=0,j,c
  
  While i<dots
    j=(i+dots-1)%dots
    ;Debug Str(i)+" "+Str(j)
    If ( ( (poly(i)\y <= *point\y) And (*point\y < poly(j)\y) ) Or
         ( (poly(j)\y <= *point\y) And (*point\y < poly(i)\y) ) ) And
       ( *point\x < (poly(j)\x-poly(i)\x) * (*point\y-poly(i)\y) / (poly(j)\y-poly(i)\y) + poly(i)\x)
      c!1
    EndIf
    i+1
  Wend
  
  ProcedureReturn c
  
EndProcedure

Procedure PointInPolygon(*polygon.TPoint, sides, *mousePosition.TPoint)
  result = #False
  
  If StartDrawing(CanvasOutput(1))
    
    result=insidePolygon(*mousePosition,polygonPoints,polygonPoints())
    
    For i = 0 To polygonPoints-1
      Plot(polygonPoints(i)\x, polygonPoints(i)\y, RGB(0,0,0))
    Next i
    
    If result
      FillArea(WindowWidth(0)/2, WindowHeight(0)/2, RGB(0,0,0), RGB(255,0,0))
    Else
      FillArea(WindowWidth(0)/2, WindowHeight(0)/2, RGB(255,255,255), RGB(255,255,255))
    EndIf
    
    DrawText(0, 0, Str(*mousePosition\x)+" | "+Str(*mousePosition\y))
    
    StopDrawing()
    Delay(1)
  EndIf
  
  ProcedureReturn result
  
EndProcedure


polygonPoints(polygonPoints)=polygonPoints(0)


wFlags = #PB_Window_SystemMenu | #PB_Window_ScreenCentered
If OpenWindow(0, 0, 0, 600, 630, "PointInPolygon()", wFlags)
  TextGadget(0, 0, 0, 600, 30, "", #PB_Text_Center)
  CanvasGadget(1, 0, 30, 600, 600, #PB_Canvas_Container)
  ;Kreis erstellen
  EinfacherKreis(polygonPoints(0), polygonPoints(1), polygonPoints(2), @Kreis)
  ;Debug "Radius: " + Kreis\r
  
  For i = 0 To 359
    polygonPoints(i)\x = ((Kreis\r) * Sin(i)) + Kreis\p\x
    polygonPoints(i)\y = ((Kreis\r) * Cos(i)) + Kreis\p\y
  Next
  
  Repeat
    event = WaitWindowEvent()
    Select event
      Case #PB_Event_CloseWindow
        appQuit = 1
      Case #PB_Event_Gadget
        If EventGadget() = 1 And EventType() = #PB_EventType_MouseMove
          mousePosition\x = GetGadgetAttribute(1, #PB_Canvas_MouseX)
          mousePosition\y = GetGadgetAttribute(1, #PB_Canvas_MouseY)
          If mouseState <> PointInPolygon(@polygonPoints(), polygonPoints, @mousePosition)
            mouseState ! 1
            If mousestate
              txt$ = "Mouse is in..."
            Else
              txt$ = "Mouse is out..."
            EndIf
            SetGadgetText(0, StrD(mousePosition\x)+" | "+StrD(mousePosition\y)+" >> "+txt$)
          EndIf
        EndIf
    EndSelect
  Until appQuit
EndIf

[chi]

Unter Windows könnte man auch CreatePolygonRgn und PtInRegion nutzen...

Code: Alles auswählen

Dim poly.POINT(2)

poly(0)\x = 10
poly(0)\y = 10
poly(1)\x = 100
poly(1)\y = 10
poly(2)\x = 50
poly(2)\y = 50

pRgn = CreatePolygonRgn_(@poly(), 3, #WINDING)

Debug PtInRegion_(pRgn, 50, 25)
Debug PtInRegion_(pRgn, 5, 5)

DeleteObject_(pRgn)