Sunday, 29 December 2013

Nested Conditional Operators

From my last post....consider this...
>>> a = [1,2,3,4,5]
>>> for i in a: True if i < 3 else ("blah" if i == 4 else False)

Which means you can extend the added function in Python 2.5 to include nested operators which is a handy function if you only have a few values to reclassify

So consider an integer field containing random numbers and an output text field, that you want to reclassify as

"True" if !Random! <  5 else "False"

would result in all the values in the Random field being reclassified to as either True or False.  If we want to expand this reclassification to include more breaks, consider the following.

"True" if !Random! <  5 else ("sort of" if !Random! == 5 else "False")

which would result in values of 5 in Random being classified as "sort of" and values > 5 as "False"

Saturday, 31 August 2013

Circular mean for directional data

When working with angles or time (minutes, days, months or years) the mean (average) of a list of values can't be determined in the same fashion as we normally do.  This example show one how to calculate the values of angles (azimuths specifically).  The function can be used in any script once the script is imported.  I will get around to other descriptive measures for circular data.


Dan Patterson
 Date  Aug 2013
  To calculate the mean (average) for points along a circle using the
  angular measure to these points from the center of the circle

  If the angles are given in degrees then there is no need to convert
  If the units are time, such as hours (24 hour clock) or months (12) then
  these need to be converted by
    angle in degrees = time quantity / number of equal intervals
    ie 0.5 = (12 hours / 24 hours per day)
    or 180 degrees
required modules
arcpy, math

def circularMean(angleList):
  angles to be given in degrees
  math module needs to be imported in the main module
  cosList = []
  cosSum = 0.0
  sinSum = 0.0
  for i in angleList:
    theCos = math.cos(math.radians(float(i)))
    theSin = math.sin(math.radians(float(i)))
    cosSum += theCos
    sinSum += theSin
  N = len(angleList)
  C = cosSum/N
  S = sinSum/N
  theMean = math.atan2(S,C)
  if theMean < 0.0:
    theMean += math.radians(360.0)
  return theMean

if __name__ == "__main__":
  import math
  sampleData = [ [350,10],
                 [350, 5],
                 [355, 10],
  print "\nCircular Mean Demo"
  print  "%20s %20s %20s" % ("angle List", "Mean radians", "Mean degrees")

  for angleList in sampleData:
    theMean = circularMean(angleList)
    print  "%20s %20.8f %20.8f" % (angleList, theMean, math.degrees(theMean))

Saturday, 24 August 2013

Segment Intersection Demo

An addendum to my previous post in the quest to find a solution to polygon diameter, antipodal distance and/or polygon fetch for concave polygons.  This ditty checks to see if line segments intersect...may come in the header etc.


  Dan Patterson
  Determines whether two line segments intersect and returns it if it does
  Detects parallel lines.  This differs from line intersections, whereby the
  point of intersection need not be on either line
Background: \
def SegmentIntersection(A,B,C,D):
  #[read the point objects, A,B,C,D which are lists of X,Y coordinates
  x0,y0 = A; x1, y1 = B; x2,y2 = C; x3,y3 = D
  # calculate differences
  Ax = x1-x0;  Ay = y1-y0
  Bx = x3-x2;  By = y3-y2
  Cx = x0-x2;  Cy = y0-y2

  # calculate the lengths of the two lines
  len1 = math.sqrt(Ax*Ax + Ay*Ay)
  len2 = math.sqrt(Bx*Bx + By*By)

  # calculate angle between the two lines
  dot = (Ax*Bx + Ay*By) #dot product
  deg = dot/(len1*len2)

  #if abs(angle) == 1:  #lines are parallel, no intersection is possible
  if(abs(deg) == 1):
    return None

  #find the intersection point (pt) between two lines
  pt = [0,0]
  div = By*Ax - Bx*Ay
  ua = (Bx*Cy - By*Cx)/div
  ub = (Ax*Cy - Ay*Cx)/div
  pt[0] = x0 + ua*Ax
  pt[1] = y0 + ua*Ay

  # calculate the combined length of the two segments
  # between pt- pnt0 and pt-pnt1
  Ax = pt[0] - x0
  Bx = pt[0] - x1
  Ay = pt[1] - y0
  By = pt[1] - y1
  segmentLen1 = math.sqrt(Ax*Ax + Ay*Ay) + math.sqrt(Bx*Bx + By*By)
  #print "segmentLen1 ", segmentLen1
  # calculate the combined length of the two segments
  # between Pt-p2 and Pt-p3
  Ax = pt[0] - x2
  Bx = pt[0] - x3
  Ay = pt[1] - y2
  By = pt[1] - y3
  segmentLen2 = math.sqrt(Ax*Ax + Ay*Ay) + math.sqrt(Bx*Bx + By*By)

  if the lengths of both sets of segments are the same as
     the lengths of the two lines the point is actually
   on the line segment.
   if the point isn't on the line, return None, else the intersection
  if(abs(len1-segmentLen1) > 0.0 or abs(len2-segmentLen2) > 0.0):
    return None
    return pt

if __name__ == "__main__":
  import math
  print "\nSegment Intersection Demo...Tests to see if two line segments intersect"
  dataList = [[[0.0,0.0],[1.0,1.0], [0.0,1.0],[1.0,0.0]],
              [[0.0,0.0], [0.0,1.0],[0.0,1.0],[1.0,1.0]],
              [[0.0,0.0], [0.0,1.0],[1.0,1.0],[0.0,1.0]],
  print dataList[0]
  for pnts in dataList:
    A,B,C,D = pnts
    pt = SegmentIntersection(A,B,C,D)
    print "\nintersection point = ", pt, " for pnts: ",A,B,C,D

Monday, 19 August 2013

Antipodal Distance Demo

This is a variant of the antipodal distance which allows you to isolate points pairs within a convex polygon forming the longest path from one side to the other side of a polygon without following the polygon edge.  If this algorithm were applied to concave polygons, then there would be no guarantee that the formed vector would be fully contained within the polygon.

Other variants of the antipodal distance can consider travel along the polygon edges, and there are many examples on the web.  Enjoy
  Dan Patterson
  To calculate the largest distance between the points the form
  a convex hull.
  - A list of list of xy pairs eg [ [1,1],[3,4,[2,6], [1,1] ]
    first and last points must be the same to close the polygon
  - numpy
  Does not account for donuts or multi-part polygons
def antipodalDistance(pnts):
  Calculates the largest distance between the points
  that make up a convex hull
  Requires: a list of list of xy pairs eg [[1,1],[3,4,[2,6],[1,1]]

  N = len(pnts)
  #thePairs = itertools.combinations(range(N), 2)   #pairs of numbers in the range N
  distList = []
  edgeList = []
  dmax = 0.0
  for i in range(0,N-2):
    for j in range(i+1, N-1):
      dist = numpy.math.hypot(pnts[i][0] - pnts[j][0], pnts[i][1] - pnts[j][1])
      diff = abs(i - j)
      if (diff == 1) or diff == abs(N-2) :
        edgeList.append( [dist, [i,j]] )
        distList.append( [dist, [i,j]] )   # compute the lower triangle of the matrix
        #distList.append( [dist, [j,i]] )   # form the upper triangle of the matrix
  dmax = distList[0]
  return [dmax, distList, edgeList]
if __name__ == "__main__":
  import numpy
  pnts = [ [0,0], [0,1], [1,4],[3,5],[5,4],[4,1],[0,0] ]  # a closed-loop convex polygon
  #pnts = [ [0,0], [2,1], [1,4],[3,5],[5,4],[4,1],[0,0] ]  # a closed-loop concave polygon

  dmax, distList, edgeList = antipodalDistance(pnts)
  print "\ndmax and index pair ", dmax
  print "\nDistList"
  for i in distList:
    print i
  print "\nedgeList:"
  for i in edgeList:
    print i

Friday, 16 August 2013

Return Azimuth...a demo script showing how to use functions

This little script shows you how to use functions to determine values and if the script is run in stand-alone mode, how default values are used.  Note the structure and requirements in the latter part of the script since this dictates the needs of the function.  There is absolutely no error checking given here.  In such cases, one would have a much longer script checking for the required inputs, to ensure that a center point is provided, a list of point coordinates are provided and to ensure that any unsavory entries (like text) are accounted for.  And for you mathphobics, study this numeric ditty.  Later

Returns an angle between 0-360
for points centered about the [0,0] origin
import math

def returnAzimuth(pntcent, pntList):
  theReturned = []
  for i in range(len(pntList)-1):
    toPnt = pntList[i]
    dY = toPnt[1] - pntCent[1]
    dX = toPnt[0] - pntCent[0]
    angle =  math.degrees(math.atan2(dY,dX))
    azimuth = 90.0 - angle                  #the difference between the x-axis and N
    azimuth = (azimuth + 360) % 360         #to account for southern/western values
    theReturned.append([pntCent, toPnt, angle, azimuth])
  return theReturned

if __name__ == "__main__":
  pntCent = [0,0]
  pntList=[ [0,1],[1,1],[1,0],[1,-1],[0,-1],[-1,-1],[-1,0],[-1,1], [-1,1] ]
  theReturned = returnAzimuth(pntCent, pntList)
  print ("\nfrom     to     angle   azimuth")
  for i in theReturned:
    print i

Tuesday, 13 August 2013

Sets continued...

I will let you figure out this one.  You should note that sets are unordered collections of data and if you want to sort an output you need to convert the set back to a list.  Here are some example of set operations.

>>> a = set( ["a","b","c"] )
>>> b = set( ["b","c","d"] )  #convert two lists to sets
>>> print a.difference(b)
>>> print b.difference(a)
>>> print a.union(b)
set(['a', 'c', 'b', 'd'])
>>> print a.intersection(b)
set(['c', 'b'])
>>> print a.symmetric_difference(b)
set(['a', 'd'])
>>> print a
set(['a', 'c', 'b'])
>>> c = list(a)
>>> c.sort()
>>> c
['a', 'b', 'c']