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)

...

True

True

False

'blah'

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"

# ObiDanGIS's Geomatics Blog

This is Dan Patterson's Geomatics blog. ObiDan is a nickname given to me by my former students. I will post ramblings about GIS and links that I find relevant to the student population and my former student-friends. Initially, I will focus on code samples that will be useful to students in the upcoming academic year. Feel free to contact me about any of the content or for help. ObiDan (ObiDanGIS@gmail.com)

## Sunday, 29 December 2013

## Saturday, 19 October 2013

### Conditional Expressions

Came across a couple of links that I will expand upon when I get some time.

http://docs.python.org/2/whatsnew/2.5.html#pep-308-conditional-expressions

http://forums.arcgis.com/threads/92241-In-line-if-statements-in-the-field-calculator-v10.0

Will greatly simplify field calculator expression.

http://docs.python.org/2/whatsnew/2.5.html#pep-308-conditional-expressions

http://forums.arcgis.com/threads/92241-In-line-if-statements-in-the-field-calculator-v10.0

Will greatly simplify field calculator expression.

## 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 CircularMean.py script is imported. I will get around to other descriptive measures for circular data.

'''

CircularMean

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Date Aug 2013

Purpose:

To calculate the mean (average) for points along a circle using the

angular measure to these points from the center of the circle

Notes:

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],

[90,270],

[180,0],

[180,360],

[0,360],

[90,180,270]

]

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))

'''

CircularMean

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Date Aug 2013

Purpose:

To calculate the mean (average) for points along a circle using the

angular measure to these points from the center of the circle

Notes:

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],

[90,270],

[180,0],

[180,360],

[0,360],

[90,180,270]

]

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 handy...read the header etc.

'''

SegmentIntersectionDemo.py

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Purpose:

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:

http://gis.stackexchange.com/questions/69025/antipodal-distance-or-polygon \

-fetch-or-polygon-diameter-for-concave-polygons

'''

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

else:

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]],

[[0.0,0.0],[0.25,0.75],[0.0,1.0],[1.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

'''

SegmentIntersectionDemo.py

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Purpose:

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:

http://gis.stackexchange.com/questions/69025/antipodal-distance-or-polygon \

-fetch-or-polygon-diameter-for-concave-polygons

'''

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

else:

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]],

[[0.0,0.0],[0.25,0.75],[0.0,1.0],[1.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

'''

AntipodalDistDemo.py

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Purpose:

To calculate the largest distance between the points the form

a convex hull.

Requirements:

- 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

Limitations:

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]] )

else:

distList.append( [dist, [i,j]] ) # compute the lower triangle of the matrix

#distList.append( [dist, [j,i]] ) # form the upper triangle of the matrix

distList.sort()

distList.reverse()

edgeList.sort()

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

Other variants of the antipodal distance can consider travel along the polygon edges, and there are many examples on the web. Enjoy

'''

AntipodalDistDemo.py

Author

Dan Patterson

obidangis@gmail.com

Dan_Patterson@carleton.ca

Purpose:

To calculate the largest distance between the points the form

a convex hull.

Requirements:

- 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

Limitations:

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]] )

else:

distList.append( [dist, [i,j]] ) # compute the lower triangle of the matrix

#distList.append( [dist, [j,i]] ) # form the upper triangle of the matrix

distList.sort()

distList.reverse()

edgeList.sort()

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

'''ReturnAzimuth

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

'''ReturnAzimuth

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)

set(['a'])

>>> print b.difference(a)

set(['d'])

>>> 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']

>>>

>>> a = set( ["a","b","c"] )

>>> b = set( ["b","c","d"] ) #convert two lists to sets

>>> print a.difference(b)

set(['a'])

>>> print b.difference(a)

set(['d'])

>>> 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']

>>>

Subscribe to:
Posts (Atom)