Estoy tratando de hacer lo que mi título dice. El objetivo es introducir una línea, segmento en la distribución uniforme de la manera por una distancia especificada, convertir los vértices a los puntos, a continuación, dibuje las líneas perpendiculares al segmento a una distancia especificada. Estas líneas, que luego se convierte en punto en una manera similar, y devolver el máximo valor de elevación de oriente y occidente de cada punto. Todo este proceso está destinado a extraer muestras de una aproximación de banco de la altura de un río. He encontrado las herramientas y los procesos de la línea de segmentación (aunque molesto que parece que uno no puede tener acceso a la elegante barra de herramientas del Editor de procesos en python), pero estoy teniendo problemas para encontrar un script que hace que el transecto perpendicular a las líneas. Nadie ha hecho algo como esto antes? Nota: he encontrado la secuencia de comandos aquí, pero tiene demasiados fallos y dejé de tratar de solucionarlo (http://forums.arcgis.com/threads/49206-Perpendicular-transects-at-regular-intervals). Gracias!
Respuestas
¿Demasiados anuncios?
Ver también Transeptos Perpendiculares a ArcGIS 10 caja de herramientas y comandos de python por Mateus Ferreira. Crea líneas perpendiculares a intervalos regulares a lo largo de la línea, pero también tiene una opción para la obtención de la división de la distancia de un nombre de campo.
La página de enlaces a varios de los debates relacionados con soluciones parciales, aquí en SIG de Intercambio de la Pila que puede ser de interés para aquellos que la construcción o la mejora de sus propias herramientas.
Mateus' script:
# Import system modules
import arcpy
from arcpy import env
import math
arcpy.env.overwriteOutput = True
#Set environments
arcpy.env.overwriteOutput = True
arcpy.env.XYResolution = "0.00001 Meters"
arcpy.env.XYTolerance = "0.0001 Meters"
# Set local variables
env.workspace = arcpy.GetParameterAsText(0)
Lines=arcpy.GetParameterAsText(1)
SplitType=arcpy.GetParameterAsText(2)
DistanceSplit=float(arcpy.GetParameterAsText(3))
TransecLength=arcpy.GetParameterAsText(4)
TransecLength_Unit=arcpy.GetParameterAsText(5)
OutputTransect=arcpy.GetParameterAsText(6)
# Def splitline module
###START SPLIT LINE CODE IN A SAME DISTANCE### Source: http://nodedangles.wordpress.com/2011/05/01/quick-dirty-arcpy-batch-splitting-polylines-to-a-specific-length/
def splitline (inFC,FCName,alongDist):
OutDir = env.workspace
outFCName = FCName
outFC = OutDir+"/"+outFCName
def distPoint(p1, p2):
calc1 = p1.X - p2.X
calc2 = p1.Y - p2.Y
return math.sqrt((calc1**2)+(calc2**2))
def midpoint(prevpoint,nextpoint,targetDist,totalDist):
newX = prevpoint.X + ((nextpoint.X - prevpoint.X) * (targetDist/totalDist))
newY = prevpoint.Y + ((nextpoint.Y - prevpoint.Y) * (targetDist/totalDist))
return arcpy.Point(newX, newY)
def splitShape(feat,splitDist):
# Count the number of points in the current multipart feature
#
partcount = feat.partCount
partnum = 0
# Enter while loop for each part in the feature (if a singlepart feature
# this will occur only once)
#
lineArray = arcpy.Array()
while partnum < partcount:
# Print the part number
#
#print "Part " + str(partnum) + ":"
part = feat.getPart(partnum)
#print part.count
totalDist = 0
pnt = part.next()
pntcount = 0
prevpoint = None
shapelist = []
# Enter while loop for each vertex
#
while pnt:
if not (prevpoint is None):
thisDist = distPoint(prevpoint,pnt)
maxAdditionalDist = splitDist - totalDist
print thisDist, totalDist, maxAdditionalDist
if (totalDist+thisDist)> splitDist:
while(totalDist+thisDist) > splitDist:
maxAdditionalDist = splitDist - totalDist
#print thisDist, totalDist, maxAdditionalDist
newpoint = midpoint(prevpoint,pnt,maxAdditionalDist,thisDist)
lineArray.add(newpoint)
shapelist.append(lineArray)
lineArray = arcpy.Array()
lineArray.add(newpoint)
prevpoint = newpoint
thisDist = distPoint(prevpoint,pnt)
totalDist = 0
lineArray.add(pnt)
totalDist+=thisDist
else:
totalDist+=thisDist
lineArray.add(pnt)
#shapelist.append(lineArray)
else:
lineArray.add(pnt)
totalDist = 0
prevpoint = pnt
pntcount += 1
pnt = part.next()
# If pnt is null, either the part is finished or there is an
# interior ring
#
if not pnt:
pnt = part.next()
if pnt:
print "Interior Ring:"
partnum += 1
if (lineArray.count > 1):
shapelist.append(lineArray)
return shapelist
if arcpy.Exists(outFC):
arcpy.Delete_management(outFC)
arcpy.Copy_management(inFC,outFC)
#origDesc = arcpy.Describe(inFC)
#sR = origDesc.spatialReference
#revDesc = arcpy.Describe(outFC)
#revDesc.ShapeFieldName
deleterows = arcpy.UpdateCursor(outFC)
for iDRow in deleterows:
deleterows.deleteRow(iDRow)
try:
del iDRow
del deleterows
except:
pass
inputRows = arcpy.SearchCursor(inFC)
outputRows = arcpy.InsertCursor(outFC)
fields = arcpy.ListFields(inFC)
numRecords = int(arcpy.GetCount_management(inFC).getOutput(0))
OnePercentThreshold = numRecords // 100
#printit(numRecords)
iCounter = 0
iCounter2 = 0
for iInRow in inputRows:
inGeom = iInRow.shape
iCounter+=1
iCounter2+=1
if (iCounter2 > (OnePercentThreshold+0)):
#printit("Processing Record "+str(iCounter) + " of "+ str(numRecords))
iCounter2=0
if (inGeom.length > alongDist):
shapeList = splitShape(iInRow.shape,alongDist)
for itmp in shapeList:
newRow = outputRows.newRow()
for ifield in fields:
if (ifield.editable):
newRow.setValue(ifield.name,iInRow.getValue(ifield.name))
newRow.shape = itmp
outputRows.insertRow(newRow)
else:
outputRows.insertRow(iInRow)
del inputRows
del outputRows
#printit("Done!")
###END SPLIT LINE CODE IN A SAME DISTANCE###
# Create "General" file geodatabase
WorkFolder=env.workspace
General_GDB=WorkFolder+"\General.gdb"
arcpy.CreateFileGDB_management(WorkFolder, "General", "CURRENT")
env.workspace=General_GDB
#Unsplit Line
LineDissolve="LineDissolve"
arcpy.Dissolve_management (Lines, LineDissolve,"", "", "SINGLE_PART")
LineSplit="LineSplit"
#Split Line
if SplitType=="Split at approximate distance":
splitline(LineDissolve, LineSplit, DistanceSplit)
else:
arcpy.SplitLine_management (LineDissolve, LineSplit)
#Add fields to LineSplit
FieldsNames=["LineID", "Direction", "Azimuth", "X_mid", "Y_mid", "AziLine_1", "AziLine_2", "Distance"]
for fn in FieldsNames:
arcpy.AddField_management (LineSplit, fn, "DOUBLE")
#Calculate Fields
CodeBlock_Direction="""def GetAzimuthPolyline(shape):
radian = math.atan((shape.lastpoint.x - shape.firstpoint.x)/(shape.lastpoint.y - shape.firstpoint.y))
degrees = radian * 180 / math.pi
return degrees"""
CodeBlock_Azimuth="""def Azimuth(direction):
if direction < 0:
azimuth = direction + 360
return azimuth
else:
return direction"""
CodeBlock_NULLS="""def findNulls(fieldValue):
if fieldValue is None:
return 0
elif fieldValue is not None:
return fieldValue"""
arcpy.CalculateField_management (LineSplit, "LineID", "!OBJECTID!", "PYTHON_9.3")
arcpy.CalculateField_management (LineSplit, "Direction", "GetAzimuthPolyline(!Shape!)", "PYTHON_9.3", CodeBlock_Direction)
arcpy.CalculateField_management (LineSplit, "Direction", "findNulls(!Direction!)", "PYTHON_9.3", CodeBlock_NULLS)
arcpy.CalculateField_management (LineSplit, "Azimuth", "Azimuth(!Direction!)", "PYTHON_9.3", CodeBlock_Azimuth)
arcpy.CalculateField_management (LineSplit, "X_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.X", "PYTHON_9.3")
arcpy.CalculateField_management (LineSplit, "Y_mid", "!Shape!.positionAlongLine(0.5,True).firstPoint.Y", "PYTHON_9.3")
CodeBlock_AziLine1="""def Azline1(azimuth):
az1 = azimuth + 90
if az1 > 360:
az1-=360
return az1
else:
return az1"""
CodeBlock_AziLine2="""def Azline2(azimuth):
az2 = azimuth - 90
if az2 < 0:
az2+=360
return az2
else:
return az2"""
arcpy.CalculateField_management (LineSplit, "AziLine_1", "Azline1(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine1)
arcpy.CalculateField_management (LineSplit, "AziLine_2", "Azline2(!Azimuth!)", "PYTHON_9.3", CodeBlock_AziLine2)
arcpy.CalculateField_management (LineSplit, "Distance", TransecLength, "PYTHON_9.3")
#Generate Azline1 and Azline2
spatial_reference=arcpy.Describe(Lines).spatialReference
Azline1="Azline1"
Azline2="Azline2"
arcpy.BearingDistanceToLine_management (LineSplit, Azline1, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_1", "DEGREES", "GEODESIC", "LineID", spatial_reference)
arcpy.BearingDistanceToLine_management (LineSplit, Azline2, "X_mid", "Y_mid", "Distance", TransecLength_Unit, "AziLine_2", "DEGREES", "GEODESIC", "LineID", spatial_reference)
#Create Azline and append Azline1 and Azline2
Azline="Azline"
arcpy.CreateFeatureclass_management(env.workspace, "Azline", "POLYLINE", "", "", "", spatial_reference)
arcpy.AddField_management (Azline, "LineID", "DOUBLE")
arcpy.Append_management ([Azline1, Azline2], Azline, "NO_TEST")
#Dissolve Azline
Azline_Dissolve="Azline_Dissolve"
arcpy.Dissolve_management (Azline, Azline_Dissolve,"LineID", "", "SINGLE_PART")
#Add Fields to Azline_Dissolve
FieldsNames2=["x_start", "y_start", "x_end", "y_end"]
for fn2 in FieldsNames2:
arcpy.AddField_management (Azline_Dissolve, fn2, "DOUBLE")
#Calculate Azline_Dissolve fields
arcpy.CalculateField_management (Azline_Dissolve, "x_start", "!Shape!.positionAlongLine(0,True).firstPoint.X", "PYTHON_9.3")
arcpy.CalculateField_management (Azline_Dissolve, "y_start", "!Shape!.positionAlongLine(0,True).firstPoint.Y", "PYTHON_9.3")
arcpy.CalculateField_management (Azline_Dissolve, "x_end", "!Shape!.positionAlongLine(1,True).firstPoint.X", "PYTHON_9.3")
arcpy.CalculateField_management (Azline_Dissolve, "y_end", "!Shape!.positionAlongLine(1,True).firstPoint.Y", "PYTHON_9.3")
#Generate output file
arcpy.XYToLine_management (Azline_Dissolve, OutputTransect,"x_start", "y_start", "x_end","y_end", "", "", spatial_reference)
#Delete General.gdb
arcpy.Delete_management(General_GDB)
He encontrado el código de la cortesía de Gerry Gabrisch publicado en un viejo ArcGIS 9.2 foro (http://arcscripts.esri.com/details.asp?dbid=15756) y modificado para que funcione en ArcGIS 10. La salida original de un archivo de texto con las coordenadas, y esta versión escribe la geometría a un vacío de shapefile. Asumir "en línea" es la entrada de la línea característica. Todo lo que necesita ser hecho una vez que se ejecuta este código es para llamar a la función "buildTransects". El código funciona bastante bien, pero todas las líneas perpendiculares en la salida no tienen la misma longitud. Tal vez necesita un poco más de ajuste, pero las longitudes de partido sólo lo suficiente para al menos discriminatoria usuario =D
# Create empty polyline for transects
spatialref = arcpy.Describe(inLine).spatialReference
if arcpy.Exists('transects'):
arcpy.Delete_management('transects')
transects = arcpy.CreateFeatureclass_management(temp,'transects.shp',"POLYLINE",'','','',spatialref)
arcpy.AddField_management(transects, "slope", "DOUBLE")
# inLine = input line, transects = new shapefile created above,
# transDist = length of transect (approximate), workspace = workspace
def buildTransects(inLine,transects,transDist,workspace):
# Define orientation (start, mid, end) and draw transect function
# This function is defined before it is called
def orientTransect(feat,ix,iy):
# If the line is horizontal or vertical, the slope and
# negative reciprocal calculations will fail, so do this instead
if starty==endy or startx==endx:
if starty == endy:
y1 = iy + transDist
y2 = iy - transDist
x1 = ix
x2 = ix
if startx == endx:
y1 = iy
y2 = iy
x1 = ix + transDist
x2 = ix - transDist
else:
# Get slope of line
m = ((starty - endy)/(startx - endx))
# Get negative reciprocal
negativereciprocal = -1*((startx - endx)/(starty - endy))
# For all values of slope, calculate perpendicular line
# with length = transDist
if m > 0:
if m >= 1:
y1 = negativereciprocal*(transDist)+ iy
y2 = negativereciprocal*(-transDist) + iy
x1 = ix + transDist
x2 = ix - transDist
if m < 1:
y1 = iy + transDist
y2 = iy - transDist
x1 = (transDist/negativereciprocal) + ix
x2 = (-transDist/negativereciprocal)+ ix
if m < 0:
if m >= -1:
y1 = iy + transDist
y2 = iy - transDist
x1 = (transDist/negativereciprocal) + ix
x2 = (-transDist/negativereciprocal)+ ix
if m < -1:
y1 = negativereciprocal*(transDist)+ iy
y2 = negativereciprocal*(-transDist) + iy
x1 = ix + transDist
x2 = ix - transDist
point1.X = x1
point1.Y = y1
point2.X = x2
point2.Y = y2
lineArray.add(point1)
lineArray.add(point2)
del x1
del x2
del y1
del y2
# Create search cursor in inLine
rows = arcpy.SearchCursor(inLine)
# Get number of records in inLine
numRecords = int(arcpy.GetCount_management(inLine).getOutput(0))
# Create new point files and array to collect values
point1 = arcpy.Point()
point2 = arcpy.Point()
lineArray = arcpy.Array()
# Define counter
counter = 0
# Loop over rows in outLine
for row in rows:
# Create the geometry object
feat = row.Shape
# Get coordinate values as lists
firstpoint = feat.firstPoint
lastpoint = feat.lastPoint
midpoint = feat.centroid
# Get X and Y values for each point
startx = firstpoint.X
starty = firstpoint.Y
endx = lastpoint.X
endy = lastpoint.Y
midx = midpoint.X
midy = midpoint.Y
m = ((starty - endy)/(startx - endx))
# For all points besides the last one
if counter < numRecords - 1:
orientTransect(feat,startx,starty)
# For the last point
else:
orientTransect(feat,endx,endy)
#Create insert cursor
cur = arcpy.InsertCursor(transects)
#Insert new row from array
feat = cur.newRow()
feat.slope = m
feat.shape = lineArray
cur.insertRow(feat)
lineArray.removeAll()
del cur
# Increase counter by 1 and start again
counter = counter + 1
del row
del rows
printit('Added %s transects to inLine and saved to temp/transects.shp' % str(counter))
Esto lo hice hace un tiempo. Esto es lo que hice. 1. Usando el algoritmo de referencia río líneas centrales, agregó escotillas de cada uno-milla (usted puede elegir el intervalo). Hechas las escotillas de 2 millas de ancho. (usted puede elegir el ancho) 2. Convierte escotillas gráfico de líneas (barra de herramientas draw) 3. Gráfico convertido líneas en una clase de entidad de línea. 4. La ruta asignada ID y vuelva a medir la característica clase de uso de "encontrar las características a lo largo de las Rutas". 5. Las líneas resultantes fueron perpendicular al río de la línea, no es necesario el general de la llanura aluvial. Girar varias secciones transversales para ajustarse a la llanura de las tendencias. Esto fue hecho en más de una ocasión. A partir de aquí yo era capaz de extraer de las elevaciones de los DEM y, a continuación, encontrar min, max, etc.
teníamos que hacer algo similar y encontró usted puede hacer manualmente. Es tedioso, pero mediante la creación de un punto al azar en la línea y, a continuación, agregar la igualdad de intervalo de puntos de ese origen, tienen iguales el intervalo de puntos. A continuación, coloca a cada punto y dibuje una línea perpendicular desde su línea de flujo o línea de vaguada, que se tiene que hacer de forma individual en cada dirección (derecha e izquierda). Usted puede combinar esas líneas perpendiculares y el uso de ellos como los transectos, o convertir a los puntos. Hay dos conjuntos de herramientas que pueden ser útiles, uno se llama ET GeoWizards que requiere una licencia, pero tiene herramientas que hacen más o menos exactamente lo que usted necesita. La otra opción es una herramienta gratuita conjunto, llamado "Geoespaciales Entorno de Modelización" (que solía ser Hawths' herramientas). El vector de herramientas de "generar puntos al azar" y "convertir las líneas de puntos". Espero que esto ayude...es una tarea tediosa.
