Start a new topic

Lookup by coordinate

Given a coordinate, what is the best way to get the polygon record at that coordinate?


Thanks! This code is helpful. I think we can get to a solution from here.


Luke

Here is a script that shows you how to do what you need. Again, I wrote this to run (and test) in Creator. You select a vertex in the scene (presumably over some terrain polygons). The script will find the polygon in the scene that is "under" this vertex. It drops the vertex to be "on the plane" of the polygon it found. It also prints out the SMC of the polygon it found.

Let us know if you have any further questions.
class NodeFound:
	def __init__(self,c):
		self.faces = list()
		self.coord = c	

def CalcPolyPlane(face):
	ok, i,j,k = mgGetPolyNormal(face)
	normal = mgvectord()
	normal.i = i
	normal.j = j
	normal.k = k
	vtx = mgGetChild(face)
	ok, x,y,z = mgGetVtxCoord(vtx)
	c = mgMakeCoord3d(x,y,z)
	plane = mgMakePlaned(c,normal)
	return plane
	
def CalcZIntersection(c, face):
	plane = CalcPolyPlane(face)
	# ax + by + cz + d = 0 
	# a,b,c,d are plane coefficients
	# x,y are x,y of coordinate
	z = (-plane.d - plane.a*c.x - plane.b*c.y)/plane.c
	return z

def GetSMC(face):
	n,c,smc = mgGetAttList(face, fltPolySmc)
	return smc

def CoordInPoly(vertx, verty, testx, testy):
# https://stackoverflow.com/questions/217578/how-can-i-determine-whether-a-2d-point-is-within-a-polygon
	nverts = len(vertx)
	c = False
	j = nverts-1
	for i in range(0,nverts):
		if ( ((verty[i]>testy) != (verty[j]>testy)) and 
				(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) ):
			c = not c
		j = i
	return c

def CoordInFaceIgnoreZ(c, face):
	vertx = list()
	verty = list()
	vtx = mgGetChild(face)
	while (vtx):
		b,x,y,z = mgGetVtxCoord(vtx)
		vertx.append(x)
		verty.append(y)
		vtx = mgGetNext(vtx)
	return CoordInPoly(vertx, verty, c.x, c.y)

def CoordInBoxIgnoreZ(c, box):
	return (
		c.x <= box.max.x and
		c.y <= box.max.y and
		c.x >= box.min.x and
		c.y >= box.min.y
		)
		
def WalkFunc (db, parent, rec, i):
	if (mgGetCode(rec) == fltPolygon):
		b, box = mgGetBounds(rec)
		if (CoordInBoxIgnoreZ(i.coord, box)):
			if (CoordInFaceIgnoreZ(i.coord, rec)):
				i.faces.append(rec)
	return MG_TRUE

db = mgGetCurrentDb()
selectList = mgGetSelectList (db)
rec,m = mgGetNextRecInList (selectList)
if (rec is None):
	mgSendMessage (MMSG_ERROR, "Nothing selected")
else:
	if (mgGetCode(rec) == fltVertex):
		b,x,y,z = mgGetVtxCoord(rec)
		c = mgMakeCoord3d(x,y,z)
		i = NodeFound(c)
		mgWalk (db, WalkFunc, None, i, MWALK_ON)
		for face in i.faces:
			print "Face:",mgGetName(face)
			mgSelectOne(face)
			z = CalcZIntersection(c,face)
			smc = GetSMC(face)
			mgSetVtxCoord(rec, c.x, c.y, z)
			print "Z intercept:",z,"SMC:",smc

 

Thanks for further describing your situation. This extra info helps and makes it a much simpler problem to solve. Presuming your terrain polygons are (generally) in the XY plane you could still use the bounding box of each polygon but when you do your comparison for "inside", ignore the Z component as all that matters is if the XY is inside the XY extents of the terrain polygon. This is like "projecting" the triangles onto the XY plane before doing the comparison. This "inside" check is still not definitive - it only gets you close to the polygon you're looking for. I assume your terrain polygons are irregular triangles so you'll still have to devise an "inside" test for whether or not the coordinate is "inside" the triangle or outside.


I hope this makes sense. I will work on another script and post it soon.

This solution may work, but perhaps if I give a broader description of the problem, then you could provide an alternate solution.


We have a set of FLT files that represent terrain maps. The terrain is stored in the FLT files as polygons. Our software is generating a set of (X,Y) coordinates (possibly millions), and we had hoped to retrieve terrain information including the Z component and SMC/FID from the FLT files at those coordinates.


Would it be feasible to take a different approach by retrieving each polygon record and computing a set of coordinates that lie on the face?


Regards,

Luke

Here is a sample script that might help you get started. I wrote this to run in Creator but if you don't have Creator you can still get the gist. When using bounding boxes on polygons, I found the mgBoxContainsCoord3d function a bit too precise for polygons that lie in any of the coordinate axes (the bounding box for such a polygons lacks 3 non-zero dimensions). I adjusted by making a very small bounding box for the coordinate (you can adjust the SmallDelta value to refine your results) and using mgBoxIntersectsBox instead.
Let us know if you have any further questions.

SmallDelta = 0.0001

def MakeCoordinateBox(c):
	box = mgboxd()
	box.min.x = c.x - SmallDelta
	box.min.y = c.y - SmallDelta
	box.min.z = c.z - SmallDelta
	box.max.x = c.x + SmallDelta
	box.max.y = c.y + SmallDelta
	box.max.z = c.z + SmallDelta
	return box

class NodeFound:
	def __init__(self,c):
		self.face = MG_NULL
		self.box = MakeCoordinateBox(c)

def WalkFunc (db, parent, rec, i):
	if (mgGetCode(rec) == fltPolygon):
		b, box = mgGetBounds(rec)
		if (mgBoxIntersectsBox(box, i.box)):
			i.face = rec
	return MG_TRUE

db = mgGetCurrentDb()
selectList = mgGetSelectList (db)
rec,m = mgGetNextRecInList (selectList)
if (rec is None):
	mgSendMessage (MMSG_ERROR, "Nothing selected")
else:
	if (mgGetCode(rec) == fltVertex):
		b,x,y,z = mgGetVtxCoord(rec)
		c = mgMakeCoord3d(x,y,z)
		i = NodeFound(c)
		print i.box.min.x,i.box.min.y,i.box.min.z
		print i.box.max.x,i.box.max.y,i.box.max.z
		mgWalk (db, WalkFunc, None, i, MWALK_ON)
		print "Face:",mgGetName(i.face)
		mgSelectOne(i.face)

 

Unfortunately the OpenFlight API does not have any kind of "spatial" search capability built in. You'd have to devise your own. 

If the coordinate you have is known to coincide exactly with a vertex of the polygon, you could mgWalk the database looking for polygons. When your walk function detects a polygon, look at each of its vertices and compare their positions with that of your coordinate. When you find the vertex whose coordinate matches, it's parent is the polygon you are looking for.

Most likely, however, the coordinate you have will NOT coincide exactly with a vertex of the polygon. If that is the case, the code gets significantly more complex.
There are some functions in the OpenFlight API that can help you:

mgGetBounds - calculates the bounding box of a node.

mgBoxContainsCoord3d - returns MG_TRUE if a coordinate is inside a bounding box
So you could mgWalk the database looking for polygons. For each polygon you hit, calculate it's bounding box. Check if the coordinate is inside that box. That will get you close. You still won't know if the coordinate lies "on the polygon". You'll have to write code to see if the coordinate is "on the same plane" as the polygon. That's pretty simple. Then you'll have to write code to see if the coordinate is on the inside of the polygon or on the outside. For convex polygons that's pretty straightforward, concave a little trickier.
I hope this can get you started. Let us know if you have further questions.

Login to post a comment