Changeset 934

Timestamp:

03/09/12 00:11:07 (15 months ago)

Author:

jjr8

Message:

Completed ArcGISPoints.CreateFeatureClassWithPoints?.

Files:

• MGET/Branches/Jason/PythonPackage/src/GeoEco/SpatialAnalysis/Points.py (modified) (4 diffs)

MGET/Branches/Jason/PythonPackage/src/GeoEco/SpatialAnalysis/Points.py

@@ -18 +18 @@
 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
+import datetime
+import math
 import os.path
 
@@ -339 +341 @@
 
     @classmethod
-    def CreatePointsAlongLines(cls, inputLines, outputFeatureClass, interval, unit, trackIDField=None, orderByField=None, startTimeField=None, endTimeField=None, pointsAtLineEnds=False, distributeRemainder=False, where=None, fieldsToCopy=None, overwriteExisting=False):
+    def CreatePointsAlongLines(cls, inputLines, outputFeatureClass, interval, unit, trackIDField=None, orderByField=None, startTimeField=None, endTimeField=None, pointsAtLineEnds=False, distributeRemainder=True, where=None, fieldsToCopy=None, overwriteExisting=False):
         cls.__doc__.Obj.ValidateMethodInvocation()
         try:
-            pass
+
+            # Perform additional validation.
+
+            isTemporal = unit in [u'days', u'hours', u'minutes', u'seconds']
+            if isTemporal:
+                if startTimeField is None:
+                    raise ValueError(_(u'A value must be provided for the Start Time Field parameter when the Interval Unit is %(unit)s.') % {u'unit': unit})
+                if endTimeField is None:
+                    raise ValueError(_(u'A value must be provided for the End Time Field parameter when the Interval Unit is %(unit)s.') % {u'unit': unit})
+
+                if trackIDField is not None and trackIDField.lower() == u'pointtime':
+                    raise ValueError(_(u'If a temporal Interval Unit is specified (i.e. days, hours, minutes, or seconds) the Track ID Field may not be named PointTime.'))
+                if orderByField is not None and orderByField.lower() == u'pointtime':
+                    raise ValueError(_(u'If a temporal Interval Unit is specified (i.e. days, hours, minutes, or seconds) the Order By Field may not be named PointTime.'))
+                if startTimeField.lower() == u'pointtime':
+                    raise ValueError(_(u'If a temporal Interval Unit is specified (i.e. days, hours, minutes, or seconds) the Start Time Field may not be named PointTime.'))
+                if endTimeField.lower() == u'pointtime':
+                    raise ValueError(_(u'If a temporal Interval Unit is specified (i.e. days, hours, minutes, or seconds) the End Time Field may not be named PointTime.'))
+                if fieldsToCopy is not None and u'pointtime' in [f.lower() for f in fieldsToCopy]:
+                    raise ValueError(_(u'If a temporal Interval Unit is specified (i.e. days, hours, minutes, or seconds) the Fields To Copy parameter may not contain a field named PointTime.'))
+
+            from GeoEco.Datasets.ArcGIS import ArcGISWorkspace, ArcGISTable
+            
+            linesTable = ArcGISTable(inputLines)
+            linesSR = linesTable.GetSpatialReference('obj')
+
+            if linesSR.IsProjected() and unit == u'decimal degrees':
+                raise ValueError(_(u'The input lines are in a projected coordinate system but the Interval Unit is decimal degrees. This is not allowed. When the lines are projected, the Interval Unit must be a linear or temporal unit (e.g. meters or days).'))
+
+            # If this the interval unit is temporal, convert the
+            # interval to a datetime.timedelta instance.
+
+            if isTemporal:
+                if unit == u'days':
+                    interval = datetime.timedelta(interval)
+                elif unit == u'hours':
+                    interval = datetime.timedelta(seconds=interval * 3600)
+                elif unit == u'minutes':
+                    interval = datetime.timedelta(seconds=interval * 60)
+                else:
+                    interval = datetime.timedelta(seconds=interval)
+
+            # Otherwise, if it is linear, convert it to meters.
+
+            elif unit != u'decimal degrees':
+                if unit == u'feet':
+                    interval *= 0.3048
+                elif unit == u'kilometers':
+                    interval *= 1000
+                elif unit == u'meters':
+                    pass
+                elif unit == u'miles':
+                    interval *= 1609.344
+                elif unit == u'nautical miles':
+                    interval *= 1852
+                elif unit == u'yards':
+                    interval *= 0.9144
+                else:
+                    raise ValueError(_(u'"%(unit)s" is an unknown linear unit. Please use feet, kilometers, meters, miles, nautical miles, or yards.') % {u'unit': unit})
+
+            # Create the output feature class.
+
+            from GeoEco.ArcGIS import GeoprocessorManager
+            from GeoEco.Datasets import QueryableAttribute
+
+            workspace = ArcGISWorkspace(os.path.dirname(outputFeatureClass), ArcGISTable, pathParsingExpressions=[u'(?P<TableName>.+)'], queryableAttributes=(QueryableAttribute(u'TableName', _(u'Table name'), UnicodeStringTypeMetadata()),))
+            outputTable = workspace.CreateTable(os.path.basename(outputFeatureClass), u'Point', linesSR)
+
+            sourceFields = [trackIDField, orderByField]
+            if fieldsToCopy is not None:
+                sourceFields.extend(fieldsToCopy)
+
+            gp = GeoprocessorManager.GetWrappedGeoprocessor()
+            destFields = []
+            i = 0
+            while i < len(sourceFields):
+                if sourceFields[i] is None:
+                    del sourceFields[i]
+                    continue
+
+                field = linesTable.GetFieldByName(sourceFields[i])
+                dataType = field.DataType
+                if dataType == u'oid':
+                    dataType = u'int32'
+                    destFields.append(sourceFields[i] + '_')
+                else:
+                    destFields.append(gp.ValidateFieldName(sourceFields[i], outputFeatureClass))
+                    
+                outputTable.AddField(destFields[i], dataType, field.Length, field.Precision, field.IsNullable)
+
+                i += 1
+
+            if isTemporal:
+                outputTable.AddField(u'PointTime', u'datetime')
+
+            # If the caller requested that the remaining distance on a
+            # track should be distributed evenly among all points,
+            # query the lines and tally the length of each track.
+
+            orderLinesBy = []
+            if trackIDField is not None:
+                orderLinesBy.append(trackIDField + ' ASC')
+            if orderByField is not None:
+                orderLinesBy.append(orderByField + ' ASC')
+            else:
+                orderLinesBy.append(inputLines.OIDFieldName + ' ASC')
+            orderLinesBy = ', '.join(orderLinesBy)
+
+            if distributeRemainder:
+                Logger.Info(_(u'Querying %(lines)s and tallying the length of each track.') % {u'lines': inputLines})
+
+                trackLengths = {}
+                gotLine = False
+                trackID = None
+                selectCursor = linesTable.OpenSelectCursor(where=where, orderBy=orderLinesBy)
+                try:
+                    while selectCursor.NextRow():
+                        if trackIDField is not None:
+                            newTrackID = selectCursor.GetValue(trackIDField)
+                            if newTrackID != trackID and gotLine:
+                                trackLengths[trackID] = trackLength
+
+                        if not gotLine or trackIDField is not None and newTrackID != trackID:
+                            gotLine = True
+                            
+                            if trackIDField is not None:
+                                trackID = newTrackID
+                            trackIsValid = True
+
+                            if isTemporal:
+                                trackLength = datetime.timedelta(0)
+                            else:
+                                trackLength = 0.
+
+                        if not trackIsValid:
+                            continue
+
+                        if isTemporal:
+                            duration = selectCursor.GetValue(endTimeField) - selectCursor.GetValue(startTimeField)
+                            if duration < datetime.timedelta(0):
+                                trackIsValid = False            # The duration of a line is not allowed to be less than zero. We will report a warning about this later.
+                                continue
+                            trackLength += duration
+                        else:
+                            trackLength += cls._GetLineLength(linesSR, unit, selectCursor, isTemporal)
+                finally:
+                    del selectCursor
+
+                if gotLine:
+                    trackLengths[trackID] = trackLength
+
+                # Now, for each track, calculate the interval that
+                # accounts for the remainder.
+
+                perTrackInterval = {}
+
+                for trackID, length in trackLengths.items():
+                    if isTemporal:
+                        intervalInSeconds = interval.microseconds * 0.000001 + interval.seconds + interval.days * 86400.
+                        lengthInSeconds = length.microseconds * 0.000001 + length.seconds + length.days * 86400.
+                        if lengthInSeconds % intervalInSeconds == 0:
+                            perTrackInterval[trackID] = datetime.timedelta(seconds = lengthInSeconds / intervalInSeconds)
+                        else:
+                            perTrackInterval[trackID] = datetime.timedelta(seconds = lengthInSeconds / (lengthInSeconds // intervalInSeconds + 1))
+
+                    elif length % interval == 0:
+                        perTrackInterval[trackID] = length / interval
+                    else:
+                        perTrackInterval[trackID] = length / (length // interval + 1)
+
+                    Logger.Debug(_(u'Track %r: length = %r, interval = %r') % (trackID, length, perTrackInterval[trackID]))
+
+            # Open an insert cursor on the output feature class and a
+            # select cursor on the lines.
+
+            Logger.Info(_(u'Inserting points into %(output)s.') % {u'output': outputFeatureClass})
+
+            insertCursor = outputTable.OpenInsertCursor()
+            try:
+                gotLine = False
+                trackID = None
+                trackIsValid = False
+                selectCursor = linesTable.OpenSelectCursor(where=where, orderBy=orderLinesBy, reportProgress=False)
+                try:
+                    while selectCursor.NextRow():
+
+                        # Read field values for this line.
+
+                        fieldValues = {}
+                        for i in range(len(sourceFields)):
+                            fieldValues[destFields[i]] = selectCursor.GetValue(sourceFields[i])
+
+                        if isTemporal:
+                            startTime = selectCursor.GetValue(startTimeField)
+                            endTime = selectCursor.GetValue(endTimeField)
+                        else:
+                            startTime = None
+                            endTime = None
+
+                        # If there are multiple tracks, get the track
+                        # ID of this line. If this line is the start
+                        # of a new track, and the caller requested
+                        # that we place lines at the end of tracks,
+                        # and we have not done so already, insert the
+                        # end point.
+                        
+                        if trackIDField is not None:
+                            newTrackID = selectCursor.GetValue(trackIDField)
+                            if newTrackID != trackID and pointsAtLineEnds and trackIsValid and (isTemporal and distanceSincePoint > datetime.timedelta(0) or not isTemporal and distanceSincePoint > 0) and endX is not None:
+                                Logger.Debug('Track %r: Writing end point' % trackID)
+                                cls._InsertPointForLine(endX, endY, endTime, insertCursor, destFields, fieldValues)
+
+                        # If this line is the start of a new track,
+                        # initialize our variables.
+
+                        if not gotLine or trackIDField is not None and newTrackID != trackID:
+                            gotLine = True
+                            
+                            if trackIDField is not None:
+                                trackID = newTrackID
+                            trackIsValid = True
+                            startOfTrack = True
+                            
+                            if isTemporal:
+                                distanceSincePoint = datetime.timedelta(0)
+                            else:
+                                distanceSincePoint = 0.
+
+                            Logger.Debug('Track %r: Started' % trackID)
+
+                            # Calculate the distance to the first
+                            # point (excluding the point that is on
+                            # the starting end of the line, if the
+                            # caller requested points on the ends). If
+                            # the caller caller requested that the
+                            # remaining distance on a track should be
+                            # distributed evenly among all points, the
+                            # the distance will be specific to this
+                            # track. If there are points at the line
+                            # ends, the distance will simply be the
+                            # track-specific point interval. If not,
+                            # then the distace will be one-half of
+                            # that interval.
+                            #
+                            # If the remaining distance should not be
+                            # distributed evenly, then the distance to
+                            # the first point is simply the point
+                            # interval requested by the caller.
+            
+                            if distributeRemainder:
+                                if pointsAtLineEnds:
+                                    distanceToNextPoint = perTrackInterval[trackID]
+                                else:
+                                    distanceToNextPoint = perTrackInterval[trackID] / 2
+                            else:
+                                distanceToNextPoint = interval
+
+                        # If the remainder of the track is invalid
+                        # (because the caller specifeid a temporal
+                        # unit and a prior line on this track had a
+                        # negative duration), skip this line.
+
+                        if not trackIsValid:
+                            Logger.Debug('Track %r: Skipping line; remainder of track is not valid.' % trackID)
+                            continue
+
+                        # If the caller specified a temporal unit,
+                        # calculate the duration of this line.
+
+                        if isTemporal:
+                            duration = endTime - startTime
+                            Logger.Debug('Track %r: Duration of line = %r, distance to next point = %r.' % (trackID, duration, distanceToNextPoint))
+
+                            # If the duration is negative, flag this
+                            # track as invalid and generate no further
+                            # points along it.
+                            
+                            if duration.seconds < 0 and duration.days < 0:
+                                if trackIDField is None:
+                                    Logger.Warning(_(u'The line with %(oid)s of %(value)r has a %(f1)s that is greater than %(f2)s. This is not allowed; it would mean the line had a negative duration. No further points will be generated.') % {u'oid': linesTable.OIDFieldName, u'value': selectCursor.GetOID(), u'f1': startTimeField, u'f2': endTimeField})
+                                else:
+                                    Logger.Warning(_(u'For the track with %(trackIDField)s of %(trackID)s, the line with %(oid)s of %(value)r has a %(f1)s that is greater than %(f2)s. This is not allowed; it would mean the line had a negative duration. No further points will not be generated along this track.') % {u'trackIDField': trackIDField, u'trackID': str(trackID), u'oid': linesTable.OIDFieldName, u'value': selectCursor.GetOID(), u'f1': startTimeField, u'f2': endTimeField})
+                                trackIsValid = False
+                                continue
+
+                            # Divide the length of this line by its
+                            # duration to get the rate of travel. Then
+                            # multiply this by the temporal distance
+                            # to the next point to get the spatial
+                            # distance along this line to that point.
+
+                            rate = cls._GetLineLength(linesSR, unit, selectCursor, isTemporal) / (duration.microseconds * 0.000001 + duration.seconds + duration.days * 86400.)
+                            spatialDistanceToNextPoint = rate * (distanceToNextPoint.microseconds * 0.000001 + distanceToNextPoint.seconds + distanceToNextPoint.days * 86400.)
+
+                            Logger.Debug('Track %r: Read new line: rate = %r, spatial distance to next point = %r' % (trackID, rate, spatialDistanceToNextPoint))
+
+                        # Iterate through the line's segments.
+
+                        geometry = selectCursor.GetGeometry()
+                        isMultipart = geometry.GetGeometryName().lower() != u'linestring'
+                        partIndex = 0
+                        endX = None
+                        endY = None
+                        
+                        while not isMultipart or isMultipart and partIndex < geometry.GetGeometryCount():
+                            if not isMultipart:
+                                part = geometry
+                            else:
+                                part = geometry.GetGeometryRef(partIndex)
+
+                            for i in range(1, part.GetPointCount()):
+
+                                # If the caller wants points at the
+                                # ends of the track and we have not
+                                # inserted the first point yet, do it.
+
+                                if startOfTrack and pointsAtLineEnds:
+                                    cls._InsertPointForLine(part.GetX(i-1), part.GetY(i-1), startTime, insertCursor, destFields, fieldValues)
+                                
+                                startOfTrack = False
+
+                                currentX = part.GetX(i-1)
+                                currentY = part.GetY(i-1)
+                                endX = part.GetX(i)
+                                endY = part.GetY(i)
+
+                                while True:
+                                    distanceToEnd = cls._GetDistance(currentX, currentY, endX, endY, linesSR, unit, isTemporal)
+                                    
+                                    if not isTemporal:
+                                        if distanceToEnd < distanceToNextPoint:
+                                            distanceSincePoint += distanceToEnd
+                                            distanceToNextPoint -= distanceToEnd
+                                            break
+
+                                        pointX, pointY = cls._GetNextPointCoords(currentX, currentY, endX, endY, distanceToNextPoint, linesSR, unit, isTemporal)
+                                        cls._InsertPointForLine(pointX, pointY, None, insertCursor, destFields, fieldValues)
+
+                                        currentX = pointX
+                                        currentY = pointY
+
+                                        distanceSincePoint = 0
+                                        if distributeRemainder:
+                                            distanceToNextPoint = perTrackInterval[trackID]
+                                        else:
+                                            distanceToNextPoint = interval
+
+                                    else:
+                                        if distanceToEnd < spatialDistanceToNextPoint:
+                                            distanceSincePoint += datetime.timedelta(seconds=distanceToEnd / rate)
+                                            distanceToNextPoint -= datetime.timedelta(seconds=distanceToEnd / rate)
+                                            spatialDistanceToNextPoint -= distanceToEnd
+                                            break
+
+                                        pointX, pointY = cls._GetNextPointCoords(currentX, currentY, endX, endY, spatialDistanceToNextPoint, linesSR, unit, isTemporal)
+                                        cls._InsertPointForLine(pointX, pointY, endTime - datetime.timedelta(seconds=distanceToEnd / rate) + distanceToNextPoint, insertCursor, destFields, fieldValues)
+
+                                        currentX = pointX
+                                        currentY = pointY
+
+                                        distanceSincePoint = datetime.timedelta(0)
+                                        if distributeRemainder:
+                                            distanceToNextPoint = perTrackInterval[trackID]
+                                        else:
+                                            distanceToNextPoint = interval
+                                        spatialDistanceToNextPoint = rate * (distanceToNextPoint.microseconds * 0.000001 + distanceToNextPoint.seconds + distanceToNextPoint.days * 86400.)
+
+                            if not isMultipart:
+                                break
+                            partIndex += 1
+
+                    # If the caller requested that we place lines at
+                    # the end of tracks, and we have not done so for
+                    # the last track, insert the end point.
+                    
+                    if pointsAtLineEnds and trackIsValid and (isTemporal and distanceSincePoint > datetime.timedelta(0) or not isTemporal and distanceSincePoint > 0) and endX is not None:
+                        Logger.Debug('Track %r: Writing end point' % trackID)
+                        cls._InsertPointForLine(endX, endY, endTime, insertCursor, destFields, fieldValues)
+
+                finally:
+                    del selectCursor
+            finally:
+                del insertCursor
         except:
             Logger.LogExceptionAsError()
             raise
+
+    @classmethod
+    def _InsertPointForLine(cls, pointX, pointY, pointTime, insertCursor, destFields, fieldValues):
+        from GeoEco.Datasets import CollectibleObject
+        insertCursor.SetGeometry(CollectibleObject._ogr().CreateGeometryFromWkt('POINT(%r %r)' % (pointX, pointY)))
+        for field in destFields:
+            insertCursor.SetValue(field, fieldValues[field])
+        if pointTime is not None:
+            insertCursor.SetValue(u'PointTime', pointTime)
+        insertCursor.InsertRow()
+
+    @classmethod
+    def _GetLineLength(cls, linesSR, unit, selectCursor, isTemporal):
+        lengthInMeters = 0.
+        geometry = selectCursor.GetGeometry()
+        isMultipart = geometry.GetGeometryName().lower() != u'linestring'
+        partIndex = 0
+
+        while not isMultipart or isMultipart and partIndex < geometry.GetGeometryCount():
+            if not isMultipart:
+                part = geometry
+            else:
+                part = geometry.GetGeometryRef(partIndex)
+
+            # If it is it is a geographic coordinate system and the
+            # unit is linear, compute the length of this part using
+            # the haversine formula.
+
+            if not linesSR.IsProjected() and unit != 'decimal degrees' and not isTemporal:
+                for i in range(1, part.GetPointCount()):
+                    lengthInMeters += cls._HaversineDistance(part.GetX(i), part.GetY(i), part.GetX(i-1), part.GetY(i-1))
+
+            # Otherwise we can just sum up the Euclidean distance of
+            # the part's segements. OGR provides a convenient function
+            # for this.
+
+            else:
+                if not linesSR.IsProjected():
+                    lengthInMeters += part.Length()
+                else:
+                    lengthInMeters += part.Length() * linesSR.GetLinearUnits()
+
+            if not isMultipart:
+                break
+            partIndex += 1
+
+        return lengthInMeters
+
+    @classmethod
+    def _GetDistance(cls, x1, y1, x2, y2, linesSR, unit, isTemporal):
+        if not linesSR.IsProjected() and unit != 'decimal degrees' and not isTemporal:
+            return cls._HaversineDistance(x1, y1, x2, y2)
+        if not linesSR.IsProjected():
+            return ((x2 - x1)**2 + (y2 - y1)**2)**0.5
+        return ((x2 - x1)**2 + (y2 - y1)**2)**0.5 * linesSR.GetLinearUnits()
+
+    @classmethod
+    def _GetNextPointCoords(cls, currentX, currentY, endX, endY, distanceToNextPoint, linesSR, unit, isTemporal):
+
+        # If it is it is a geographic coordinate system and the unit
+        # is linear, we have to perform this calculation on a great
+        # circle. First calculate the forward azimuth (the initial
+        # bearing) toward the end point.
+
+        if not linesSR.IsProjected() and unit != 'decimal degrees' and not isTemporal:
+            lat1 = math.radians(currentY)
+            lat2 = math.radians(endY)
+            dlon = math.radians(endX - currentX)
+
+            y = math.sin(dlon) * math.cos(lat2)
+            x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(dlon)
+            azimuth = math.atan2(y, x)
+
+            # Now compute the coordinates of the next point, traveling
+            # the specified distance on the forward azimuth.
+
+            distRadians = distanceToNextPoint / 6371009     # 6371009 is mean radius of Earth, in meters
+
+            lon1 = math.radians(currentX)
+
+            lat2 = math.asin(math.sin(lat1) * math.cos(distRadians) + math.cos(lat1) * math.sin(distRadians) * math.cos(azimuth))
+            lon2 = lon1 + math.atan2(math.sin(azimuth) * math.sin(distRadians) * math.cos(lat1), math.cos(distRadians) - math.sin(lat1) * math.sin(lat2))
+            lon2 = (lon2 + 3 * math.pi) % (2 * math.pi) - math.pi;      # Normalize to -180 to +180
+
+            return math.degrees(lon2), math.degrees(lat2)
+
+        # Otherwise we can perform this calculation using basic
+        # trigonometry.
+
+        proportion = distanceToNextPoint / (((currentX - endX)**2 + (currentY - endY)**2)**0.5)
+        return currentX + (endX - currentX)*proportion, currentY + (endY - currentY)*proportion
+
+    @classmethod
+    def _HaversineDistance(cls, lon1, lat1, lon2, lat2):
+        lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2])
+        dlon = lon2 - lon1
+        dlat = lat2 - lat1
+        a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
+        c = 2 * math.asin(math.sqrt(a))
+        return 6371009 * c      # 6371009 is mean radius of Earth, in meters
 
     @classmethod
@@ -639 +1123 @@
             nearAreProjected = nearFeaturesDescribe.SpatialReference.ProjectionCode != 0
 
-        # Validate that either the points or the near features are
+        # Warn the user if neither the points or the near features are
         # projected.
 
         if not pointsAreProjected and not nearAreProjected:
-            Logger.RaiseException(ValueError(_(u'Neither the point features (%(points)s) nor the near features (%(near)s) use a projected coordinate system. One or the other of these must use a projected coordinate system. Please project one or the other of them and try again.') % {u'points': pointFeatures, u'near': nearFeatures}))
+            Logger.Warning(_(u'Neither the point features (%(points)s) nor the near features (%(near)s) use a projected coordinate system. Because ArcGIS finds the nearest feature in Euclidean distance, not geodetic distance, the features identified as nearest may not actually be the nearest geodetically.') % {u'points': pointFeatures, u'near': nearFeatures})
 
         # Copy the caller's points into the temp geodatabase,
@@ -1768 +2252 @@
 interval. If the lines are in a geographic coordinate system, points
 will be placed along the great circles that connect the lines'
-vertices and spaced according to geodesic distance computed on the
-coordiante system's ellipsoid. Because of this, the points will not
-occur exactly along the lines unless the lines are themselves great
-circles (i.e. they exactly follow the equator or a meridian).
+vertices and spaced according to haversine distance. Because of this,
+the points will not occur exactly along the lines unless the lines are
+themselves great circles (i.e. they exactly follow the equator or a
+meridian).
 
 If a temporal unit is used, the input lines may be in a projected or