You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
If you or someone in your organization still prefers to draw cross sections "the old fashioned way" then this script is for you!
Supply a feature class with a single polyline and specify a DEM (GRID, tiff, img, etc.; NOT compatible with mosaic datasets... yet) and it generates a PolylineZ (3D) feature class and a 2D profile plot and exports it to an SVG file. Note that the spatial references of all inputs must match and the script assumes that the x, y, and z values are in meters. You can specify the horizontal and vertical grid spacing of the plot as well as the vertical exaggeration and map scale. Optionally annotates the profile with info from the user-specified GeMS feature classes. The SVG file is written using the Inkscape SVG format which allows the use of "layers" that Illustrator interprets as groups, keeping the grid, profile, and annotations organized and easy to separate.
I had created a basic profile graphic tool many years ago to solve the problem that the ArcMap 3D/Spatial Analyst profile tool produced an accurate but unsatisfactory profile graphic. The ArcMap output was not scaled to the map layout and the horizontal and vertical scales were in no way related and none of these options could be specified which required scaling the plot later in Illustrator. Then the geologic map information had to be added manually by superimposing the profile on the map and transferring the data to the profile. My script created a scaled profile, but it did not add any geologic information as the lack of a common schema to the inputs made this difficult to automate.
After a few fits and starts, I finally finished this script when I had to produce new profiles for one of our geologists to complete some cross sections for a previously unpublished map that I was compiling. Drafts of the cross sections had been hand drawn, but the profiles were generalized a bit too much for our liking and the contacts were adjusted during the compilation process making the hand drawn cross sections obsolete.
Key features of this script:
NO dependency on any ArcGIS extensions (3D/Spatial Analyst).
It is much faster at interpolating a 3D feature than the built-in arcpy interpolation routine (arcpy.management.GetCellValue) by using a numpy array (arcpy.RasterToNumPyArray) to store a "clip" of the DEM that matches the extent of the input line and reading the Z values from that array instead of querying the raster dataset itself.
The sample interval can be set independent of raster cell size or vertex location on the line (which some might consider a problem?)
The script can generate a profile to match any map scale at the desired vertical exaggeration.
Check out the code below for a detailed description and see an example of the output in the attached SVG file. The example is a re-creation of the profile shown in cross section B-B' on the Bacon Gap geologic quadrangle map: https://ngmdb.usgs.gov/Prodesc/proddesc_108498.htm. The cross section on the published map was created by the author; I only re-created the profile and plotted the pertinent geologic map data to show how the tool output can be used as the starting point for drawing a cross section.
I have also attached an ArcGIS Pro .atbx file with the script in a zip file to make setting up the tool easier since it has a lot of parameters. Unzip, add the .atbx file to your ArcGIS Pro Favorites, verify in the script properties that the path to the .py file is correct, and run the tool. The tool does have one additional Python dependency (svgwrite) that you may have to add to your Python environment if you don't already have it. In the folder of your cloned ArcGIS Python environment, install svgwrite (Start > Python command prompt > pip install svgwrite) and the code should work when run from an ArcGIS session.
"""
GeneratePolylineZExport2DProfileSVG.py
Andrew L. Wunderlich
andrew.wunderlich@tn.gov
2024-10-30
Takes a single line feature and a DEM and produces a z-enabled polyline and draws a 2D profile of the line on a
measured grid based on the input parameters of map scale and vertical exaggeration. Optionally, the script can add
additional annotations (symbols and labels) based on GeMS feature classes ContactsAndFaults, GeologicLines,
MapUnitLines, MapUnitPolys, and OrientationPoints.
Toolbox setup:
General:
Name: GeneratePolylineZAndExportProfileToSVGFile
Label: Generate PolylineZ and export profile to SVG file
Description: Generates a PolylineZ (3D) feature class from a polyline and a DEM and generates a 2D profile and exports
it to an SVG file. Optionally annotates the profile with info from the user-specified GeMS feature classes.
Parameters:
Label Name Data Type Type Direction Filter Dependency Default
0 Workspace output_ws Workspace Required Input
1 Input polyline feature class input_fc Feature Class Required Input Feature Type: Polyline output_ws
2 DEM dem Raster Dataset Required Input
3 Sample interval in map units sample_interval Double Required Input
4 Output ZM polyline feature class name output_fc String Required Output
5 Output SVG file output_svg File Required Output File
6 SVG grid interval (X) in map units interval_length Double Required Input
7 SVG grid interval (Z) in map units interval_elevation Double Required Input
8 Map scale denominator (e.g. 24000) map_scale Double Required Input
9 Vertical exaggeration (1 = no exaggeration) vertical_ exaggeration Double Required Input 1
10 Annotate profile? annotate_profile Boolean Optional Input FALSE
11 ContactsAndFaults ContactsAndFaults Feature Class Optional Input Feature Type: Polyline
12 GeologicLines Geologiclines Feature Class Optional Input Feature Type: Polyline
13 MapUnitLines MapUnitlines Feature Class Optional Input Feature Type: Polyline
14 MapUnitPolys MapUnitPolys Feature Class Optional Input Feature Type: Polygon
15 OrientationPoints OrientationPoints Feature Class Optional Input Feature Type: Point
16 Search distance search_distance Double Optional Input
Validation:
def updateParameters(self):
# Modify parameter values and properties.
# This gets called each time a parameter is modified, before
# standard validation.
if str(self.params[10].value).upper() == "TRUE":
self.params[11].enabled = True
self.params[12].enabled = True
self.params[13].enabled = True
self.params[14].enabled = True
self.params[15].enabled = True
self.params[16].enabled = True
else:
self.params[11].enabled = False
self.params[12].enabled = False
self.params[13].enabled = False
self.params[14].enabled = False
self.params[15].enabled = False
self.params[16].enabled = False
return
Dependencies:
-arcpy: To interact with the PolylineZM feature class and extract the geometry and Z values.
-svgwrite: To generate the SVG graphic. You can install it using pip install svgwrite.
Code explanation:
PolylineZ Portion
The basic process involves:
-Extracting points along the polyline at intervals equal to the raster DEM's cell size.
-Sampling the DEM at these points to get the Z (elevation) values.
-Constructing a new PolylineZM feature class using the extracted points and their corresponding Z values.
Script Explanation:
1. Extract Points Along the Polyline (extract_z_from_dem):
-We use positionAlongLine to extract points along the polyline at regular intervals based on the cell size of the DEM.
-arcpy.GetCellValue_management is used to get the Z (elevation) values from the DEM at the extracted point locations.
2. Create an Empty PolylineZM Feature Class (create_polyline_zm_feature_class):
-A new feature class with polyline geometry and Z-enabled fields is created using arcpy.management.CreateFeatureclass.
3. Add PolylineZM Features (add_polyline_zm_to_feature_class):
-A polyline is created using the extracted points and their corresponding Z values.
-arcpy.da.InsertCursor is used to insert this polyline into the new feature class.
4. Processing Each Polyline (process_polyline_zm):
-The script processes each polyline feature in the input feature class, extracts Z values along it, and adds the
polyline with Z values to the new PolylineZM feature class.
Input Parameters:
1. input_fc: The input polyline feature class.
2. dem: The DEM raster file.
3. output_fc: The output PolylineZM feature class.
Notes:
-This script assumes that the DEM and polyline feature class share the same spatial reference or that appropriate
transformations are handled in advance.
-This script assumes units of meters for the x, y, and z values
-It uses the user-specified value (in map units) to determine the interval at which to extract points along the
polyline.
-The script assumes that the environment is set up to handle file geodatabases.
SVG Draw/Export Portion
The basic process involves:
-Use arcpy to read the feature class geometry and Z values.
-Convert the PolylineZM's 3D data into 2D data (Z values for the y-axis, and cumulative length for the x-axis).
-Use svgwrite to generate a scalable vector graphic (SVG).
Explanation:
1. Extracting Geometry with arcpy.da.SearchCursor:
Use arcpy to extract PolylineZM vertices and their associated Z values.
-The script uses the SHAPE@ token to retrieve the geometry from the feature class. It then processes each vertex in
the polyline, computing the cumulative distance along the line (x-axis) and recording the Z value for each vertex
(y-axis).
2. Scaling and Plotting:
Plot the data using a simple 2D SVG graphic with labeled axes, where the x-axis represents the length along the
polyline and the y-axis represents elevation (Z):
-The length_scale is calculated based on the total length of the polyline and the width of the plot area.
-The elevation_scale is calculated using the range of Z values and the height of the plot area.
-The line is plotted as a polyline in the SVG, with x-coordinates representing the cumulative distance and
y-coordinates representing the elevation values, which are scaled by the user-specified vertical exaggeration.
3. Grid and Labels:
Add grid lines and labels based on user-specified intervals for distance and elevation:
-Horizontal grid lines are added at intervals of the user-specified length (x-axis).
-Vertical grid lines are added at intervals of the user-specified elevation (y-axis).
-Labels are added to both axes for reference.
4. SVG Output:
-The svgwrite library is used to generate the SVG output, creating a scalable vector graphic that can be viewed in any
SVG-compatible viewer or browser.
Customization:
-You can adjust the plot grid spacing interval_length (X) and interval_elevation (Z) based on your data and scale.
-The plot_width and plot_height values in the script are controlled by the map_scale value.
-The vertical scale is controlled by the vertical_ex value.
Annotation Portion
The basic process involves:
-Intersect the GeMS feature classes with the profile line
-Read the attributes and symbolize and label the profile plot accordingly
Explanation:
1. Intersecting line and polygon geometry:
-Intersects the geometry of each GeMS fc with the profile line using polyline.intersect which returns the topology of
the lowest order, i.e., line-line => point, line-polygon => line
-Reads selected attributes of the intersected feature and plots them based on the location along the profile
-Polygons that intersect the profile more than once result in a multipart polyline that must be split and each part's
location calculated to accurately label the profile.
-Symbols on the profile are plotted and labeled based on Type and Label attributes
2. Find and project point features to the profile line:
-Creates a temporary point feature class to store the projected intersection points
-Searches for points within the user-specified radius search_distance
-Uses Pythagorean theorem to estimate the location of the point along the line
-Plots a symbol and "rotates" it using the inclination and azimuth to indicate the direction of dip along the line
"""
import sys
import arcpy
import os
import svgwrite
from svgwrite.extensions import Inkscape
from math import sqrt, radians, cos, sin
import datetime
# Copied form GeMS Utility scripts
def addMsgAndPrint(msg, severity=0):
# prints msg to screen and adds msg to the geoprocessor (in case this is run as a tool)
# print msg
try:
for string in msg.split("\n"):
# Add appropriate geoprocessing message
if severity == 0:
arcpy.AddMessage(string)
elif severity == 1:
arcpy.AddWarning(string)
elif severity == 2:
arcpy.AddError(string)
except:
pass
### THE NEXT SIX SCRIPTS WERE COPIED AND ADAPTED FROM CODE BY XANDER BAKKER (2014)
### https://community.esri.com/t5/python-documents/extract-raster-values-using-numpy/ta-p/907274
def getListOfPointsFromPolyline(line, interval):
points = []
d = 0
max_len = line.length
while d < max_len:
point = line.positionAlongLine(d)
points.append(point.centroid)
d = d + interval
points.append(line.lastPoint)
return points
def getRowColExtentFromXY(ext, point, cell_size):
# row, col start at 0 from upper left
row = int(((ext.YMax - point.Y) - ((ext.YMax - point.Y) % cell_size)) / cell_size)
col = int(((point.X - ext.XMin) - ((point.X - ext.XMin) % cell_size)) / cell_size)
return row, col
def createNumpyArrayUsingExtent(raster, ext, cell_size):
lowerLeft = arcpy.Point(ext.XMin, ext.YMin)
# Added and extra col and row to the array to avoid out of range error on endpoint (last PolylineZ extraction point)
nCols = int(ext.width / cell_size) + 1
nRows = int(ext.height / cell_size) + 1
addMsgAndPrint(f"Array of {nRows} "
f"rows and {nCols} "
f"columns specified for extraction from DEM")
try:
return arcpy.RasterToNumPyArray(raster, lowerLeft, nCols, nRows, -9999)
except RuntimeError:
e = sys.exc_info()[1]
addMsgAndPrint(f"The array of {nRows} rows and {nCols} columns could not be extracted. " + e.args[0], 2)
exit()
def adaptExtentUsingRaster(ext, raster, cell_size):
# This tool does not work with mosaic datasets
ras_ext = arcpy.Describe(raster).extent
xMin = ext.XMin - ((ext.XMin - ras_ext.XMin) % cell_size)
yMin = ext.YMin - ((ext.YMin - ras_ext.YMin) % cell_size)
xMax = ext.XMax + ((ras_ext.XMax - ext.XMax) % cell_size)
yMax = ext.YMax + ((ras_ext.YMax - ext.YMax) % cell_size)
return arcpy.Extent(xMin, yMin, xMax, yMax)
def getRasterCellSize(raster):
desc = arcpy.Describe(raster)
return (desc.meanCellHeight + desc.meanCellWidth) / 2
def getExtentFromPolyline(line):
return line.extent
### ADAPTED FROM SCRIPT BY XANDER BAKKER AND SUPPLEMENTED BY ALW
### https://community.esri.com/t5/python-documents/extract-raster-values-using-numpy/ta-p/907274
def extract_z_from_dem(dem, polyline, interval):
# Extract Z values from DEM along a polyline at specified intervals.
# Determine raster cell size for sampling
cell_size = getRasterCellSize(dem)
addMsgAndPrint(f"Cell size of DEM (in units of projection): {str(cell_size)}")
addMsgAndPrint(f"Sampling interval (in units of projection): {str(interval)}")
# Get extent from line and adapt extent to raster
ext1 = getExtentFromPolyline(polyline)
addMsgAndPrint(f"Extent of polyline: "
f"XMin {str(ext1.XMin)}, "
f"YMin {str(ext1.YMin)}, "
f"XMax {str(ext1.XMax)}, "
f"YMax {str(ext1.YMax)}")
ext2 = adaptExtentUsingRaster(ext1, dem, cell_size)
addMsgAndPrint(f"Extent of DEM (adapted from polyline) to be extracted to array: "
f"XMin {str(ext2.XMin)}, "
f"YMin {str(ext2.YMin)}, "
f"XMax {str(ext2.XMax)}, "
f"YMax {str(ext2.YMax)}")
# Create numpy array
npArray = createNumpyArrayUsingExtent(dem, ext2, cell_size)
"""addMsgAndPrint(f"Array of {str(npArray.shape[0])} "
f"rows and {str(npArray.shape[1])} "
f"columns extracted from DEM using adapted polyline extent")"""
# Generate points along the polyline at the specified interval
# [points] is for building the new PolylineZ (Array of point x, y coordinate pairs)
points = []
# Populate [points] array with x, y coordinate pairs
current_distance = 0
while current_distance < polyline.length:
point = polyline.positionAlongLine(current_distance)
points.append((point.getPart().X, point.getPart().Y))
current_distance += interval
points.append((polyline.lastPoint.X, polyline.lastPoint.Y))
# [pointsZ] is for extracting the Z values from the DEM (Array of point geometries)
pointsZ = getListOfPointsFromPolyline(polyline, interval)
# Get the Z values from the DEM for each point and print the info
addMsgAndPrint("Data table of Z values for points sampled along the polyline:")
addMsgAndPrint("interval\traster_row\traster_col\tpoint_X\tpoint_Y\tZ_value")
z_values = []
current_distance = 0
for point in pointsZ:
try:
row, col = getRowColExtentFromXY(ext2, point, cell_size)
z_value = npArray.item(row, col)
addMsgAndPrint(f"{current_distance}\t{row}\t{col}\t{point.X}\t{point.Y}\t{z_value}")
current_distance += interval
z_values.append(z_value)
except IndexError:
addMsgAndPrint(f"The item at index {row}, {col} is out of range. Moving on...")
break
return points, z_values
def create_polyline_zm_feature_class(output_fc_name, spatial_ref):
# Create a new empty PolylineZM feature class.
arcpy.CreateFeatureclass_management(out_path=output_ws,
out_name=output_fc_name,
geometry_type="POLYLINE",
has_m="DISABLED",
has_z="ENABLED",
spatial_reference=spatial_ref)
addMsgAndPrint(f"New (empty) PolylineZM feature class {output_fc_name} created.")
def add_polyline_zm_to_feature_class(output_fc, points, z_values, spatial_ref):
#Add a polyline with Z values to a PolylineZM feature class.
# Create array to hold the points
array = arcpy.Array()
for (x, y), z in zip(points, z_values):
point = arcpy.Point(x, y, z) # Create a Point with Z
array.add(point)
# Create a PolylineZM geometry
polyline_zm = arcpy.Polyline(array, spatial_ref, has_z=True, has_m=True)
# Insert the polyline into the feature class
with arcpy.da.InsertCursor(output_fc, ["SHAPE@"]) as cursor:
cursor.insertRow([polyline_zm])
addMsgAndPrint(f"XYZ point array created for input polyline {os.path.split(input_fc)[1]}.")
def process_polyline_zm(input_fc, dem, feature_class, spatial_ref):
# Main function to extract Z values along a polyline from DEM and create a PolylineZM feature class.
# Get the cell size (resolution) of the DEM
#desc = arcpy.Describe(dem)
#cell_size = desc.meanCellHeight # Assuming square cells
# Use sample interval specified by user instead
interval = float(sample_interval)
# Create an empty PolylineZM feature class
create_polyline_zm_feature_class(output_fc_name, spatial_ref)
# Iterate through each polyline in the input feature class
with arcpy.da.SearchCursor(input_fc, ["SHAPE@"]) as cursor:
for row in cursor:
polyline = row[0]
#points, z_values = extract_z_from_dem(dem, polyline, cell_size)
# Use sample interval specified by user instead
points, z_values = extract_z_from_dem(dem, polyline, interval)
# Add the new polyline with Z values to the output feature class
add_polyline_zm_to_feature_class(feature_class, points, z_values, spatial_ref)
# Let user know process was successful
addMsgAndPrint(f"PolylineZ feature added to feature class {output_fc_name}.")
### SVG PART OF SCRIPT ###
def calculate_distance(point1, point2):
# Calculate the 2D distance between two points (ignoring Z).
return sqrt((point2.X - point1.X) ** 2 + (point2.Y - point1.Y) ** 2)
#return sqrt((point2.getPart().X - point1.getPart().X) ** 2 + (point2.getPart().Y - point1.getPart().Y) ** 2)
def plot_elevation_profile(zm_fc, output_svg, interval_length, interval_elevation, map_scale, vertical_ex, contacts_fc,
geolines_fc, mulines_fc, mupolys_fc, orpoints_fc):
# Variables for scaling the SVG plot
ppi = 72 # pixels per inch in an SVG file
meters_inch = 0.0254 # meters per inch
total_length = 0 # Total length along the polyline (assumed meters)
min_z = float('inf')
max_z = float('-inf')
# List to store the cumulative distance and elevation points
profile_points = []
# Fetch the first polyline geometry from the feature class
with arcpy.da.SearchCursor(zm_fc, ["SHAPE@"]) as cursor:
for row in cursor:
polyline = row[0] # Assuming single polyline
# For each part in the polyline
for part in polyline:
previous_point = None
cumulative_length = 0
# For each point in the part
for point in part:
if previous_point:
distance = calculate_distance(previous_point, point)
cumulative_length += distance
else:
cumulative_length = 0
# Update min and max Z values ignoring NoData value in the array (-9999)
if point.Z is not None and point.Z != -9999:
min_z = min(min_z, point.Z)
max_z = max(max_z, point.Z)
# Store cumulative length and Z value
profile_points.append((cumulative_length, point.Z))
previous_point = point
total_length = cumulative_length
# Establish the Z range for the profile plot in map units
z_range = max_z - min_z
# Establish the height of the Z grid plot (grid starts from 0; "sea level")
z_range_grid = max_z
# Determine the size of drawing based on length of profile, map scale, and vertical exaggeration
# Width of drawing
mu_inch = float(map_scale) * meters_inch
plot_width_inch = total_length / mu_inch
# Height of drawing considering vertical exaggeration
"""plot_height_inch = (z_range / float(mu_inch)) * float(vertical_ex)"""
plot_height_inch = (z_range_grid / float(mu_inch)) * float(vertical_ex)
# Define margins and dimensions for the plot area in pixels
margin = 50 # pixels
plot_width = plot_width_inch * ppi # convert inches to pixels
plot_height = plot_height_inch * ppi # convert inches to pixels
# Calculate scaling factors for the SVG plot
length_scale = plot_width / total_length
"""elevation_scale = plot_height / z_range if z_range > 0 else 1"""
elevation_scale = plot_height / z_range_grid if z_range > 0 else 1
# Added for new plotting method below
profile_y_base = margin + plot_height
# Create an SVG drawing
"""dwg = svgwrite.Drawing(output_svg, size=("1000px", "500px"))"""
dwg_width = plot_width + (margin * 2)
dwg_height = plot_height + (margin * 2)
dwg = svgwrite.Drawing(output_svg, size=(dwg_width, dwg_height))
addMsgAndPrint(f"Creating {os.path.split(output_svg)[1]}; "
f"page size {round(dwg_width / 72, 4)} in. w x {round(dwg_height / 72, 4)} in. h.")
# Activate the Inkscape extension of the SVG drawing for building "layers" (groups in Illustrator)
inkscape = Inkscape(dwg)
# Create "layers" in the Inkscape SVG for each of the drawing elements
grid_layer = inkscape.layer(label="Grid", locked=False)
misc_layer = inkscape.layer(label="Doc Info", locked=False)
profile_layer = inkscape.layer(label="Profile", locked=False)
ticks_layer = inkscape.layer(label="CAF ticks", locked=False)
geolines_layer = inkscape.layer(label="GEL ticks", locked=False)
mulines_layer = inkscape.layer(label="MUL ticks", locked=False)
mupolys_layer = inkscape.layer(label="MUP labels", locked=False)
dip_layer = inkscape.layer(label="ORP symbols", locked=False)
# Add axes and grid lines
dwg.add(grid_layer)
# X axis grid
for i in range(0, int(total_length), int(interval_length)):
# Define the plot location
x = margin + i * length_scale
# Draw the vertical line
grid_layer.add(dwg.line(start=(x, margin), end=(x, margin + plot_height), stroke="lightgray", stroke_width=0.5))
# Add the label
grid_layer.add(dwg.text(f"{i:,}", insert=(x, margin + plot_height + 12), fill="black", font_size="8px",
text_anchor="middle"))
# Y axis grid
for i in range(0, int(max_z), int(interval_elevation)):
# Define the plot location
y = (margin + plot_height) - i * elevation_scale
# Draw the horizontal line
grid_layer.add(dwg.line(start=(margin, y), end=(margin + plot_width, y), stroke="lightgray", stroke_width=0.5))
# Add the label
grid_layer.add(
dwg.text(f"{i:,}", insert=(margin - 8, y + 3.6), fill="black", font_size="8px", text_anchor="end"))
addMsgAndPrint(f" Measured grid written to SVG file...")
# Add some additional labels to the plot to identify key values
dwg.add(misc_layer)
# Title
misc_layer.add(
dwg.text(f"{os.path.split(output_svg)[-1]}", insert=(dwg_width / 2, margin / 2), fill="black", font_size="10px",
text_anchor="middle"))
# Tool parameters
if ticks_and_anno == "true":
misc_layer.add(
dwg.text(
f"Output scale 1:{int(map_scale):,}; "
f"Vertical exaggeration: {vertical_ex}x; "
f"Sample interval: {sample_interval} m; "
f"ORP search distance: {search_distance} m; "
f"DEM: {os.path.split(dem)[-1]}",
insert=(dwg_width - margin, dwg_height - (margin / 3) - 7.2),
fill="black", font_size="6px", text_anchor="end"))
else:
misc_layer.add(
dwg.text(
f"Output scale 1:{int(map_scale):,}; "
f"Vertical exaggeration: {vertical_ex}x; "
f"Sample interval: {sample_interval} m; "
f"DEM: {os.path.split(dem)[-1]}",
insert=(dwg_width - margin, dwg_height - (margin / 3) - 7.2),
fill="black", font_size="6px", text_anchor="end"))
# Z profile info
misc_layer.add(
dwg.text(f"ZMin: {round(min_z, 3):,} m, "
f"ZMax: {round(max_z, 3):,} m, "
f"Length: {round(total_length, 3):,} m",
insert=(dwg_width - margin, dwg_height - (margin / 3)),
fill="black", font_size="6px", text_anchor="end"))
addMsgAndPrint(f" Document and profile info written to SVG file...")
# Function to map profile coordinates
def map_to_svg_coords(dist, z):
x = margin + dist * length_scale
#y = profile_y_base - (z - min_z) * elevation_scale
y = profile_y_base - z * elevation_scale
#y = z * elevation_scale
return x, y
# Function to find closest Z value on profile to the dist value of a point along the profile
def find_closest_Z_value(data, target):
# Find the tuple with the closest first value to the target
closest = min(data, key=lambda x: abs(x[0] - target))
# Return the second value in the closest tuple
return closest[1]
# Draw the elevation profile
dwg.add(profile_layer)
dwgPoints = [map_to_svg_coords(dist, z) for dist, z in profile_points if z is not None]
profile_layer.add(dwg.polyline(dwgPoints, stroke="blue", fill="none", stroke_width=0.4252))
addMsgAndPrint(f" Profile line written to SVG file...")
if ticks_and_anno == "true":
addMsgAndPrint(f"Option to annotate profile from geologic data selected. Gathering info...")
# Create SVG layers to hold these additional items
dwg.add(ticks_layer)
dwg.add(geolines_layer)
dwg.add(mulines_layer)
dwg.add(mupolys_layer)
dwg.add(dip_layer)
# Process line intersections for CAF tick marks
if len(contacts_fc) > 3:
with arcpy.da.SearchCursor(contacts_fc, ["SHAPE@", "Type", "Label"]) as cursor:
for row in cursor:
line_geom, line_type, line_label = row
intersection_points = polyline.intersect(line_geom, 1)
for point in intersection_points:
"""addMsgAndPrint(f"Map coords of CAF-profile intersection: X: {point.X}, Y: {point.Y}")"""
dist = calculate_distance(polyline.firstPoint, point)
x, y = map_to_svg_coords(dist, point.Z) # Location on 2D profile
txt_rotate = "rotate(270,%s, %s)" % (x, y) # Set text rotation to vertical
# Define the tick as a polyline with these vertices
line_points = [(x, y + 3), (x, y), (x, y - 3)]
# Symbolize different line types (contacts, faults, etc.)
if str(line_type)[:5].lower() == "fault":
ticks_layer.add(dwg.polyline(line_points, stroke="red", stroke_width=1.0630))
# Label the tick
if line_label:
ticks_layer.add(dwg.text(line_label, insert=(x + 10, y + 1.5),
font_size="6px", fill="black", transform=txt_rotate))
elif str(line_type)[:7].lower() == "contact":
ticks_layer.add(dwg.polyline(line_points, stroke="black", stroke_width=0.4252))
# Label the tick
"""if len(line_label) > 0:
ticks_layer.add(dwg.text(line_label, insert=(x + 10, y + 1.5),
font_size="6px", fill="black", transform=txt_rotate))"""
else:
ticks_layer.add(dwg.polyline(line_points, stroke="yellow", stroke_width=0.5))
# Only label the tick if it is not a contact or fault
if line_label:
ticks_layer.add(dwg.text(line_type + ", " + line_label, insert=(x + 5, y + 1.5),
font_size="6px", fill="black", transform=txt_rotate))
else:
ticks_layer.add(dwg.text(line_type, insert=(x + 5, y + 1.5),
font_size="6px", fill="black", transform=txt_rotate))
addMsgAndPrint(f" CAF intersections written to SVG file...")
# Process line intersections for GEL tick marks
if len(geolines_fc) > 3:
with arcpy.da.SearchCursor(geolines_fc, ["SHAPE@", "Type", "Label"]) as cursor:
for row in cursor:
line_geom, line_type, line_label = row
intersection_points = polyline.intersect(line_geom, 1)
for point in intersection_points:
"""addMsgAndPrint(f"Map coords of GEL-profile intersection: X: {point.X}, Y: {point.Y}")"""
dist = calculate_distance(polyline.firstPoint, point)
x, y = map_to_svg_coords(dist, point.Z) # Location on 2D profile
# Define the tick as a polyline with these vertices
line_points = [(x, y + 3), (x, y), (x, y - 3)]
# Define rotations for fold symbols and labels
anti_rot = "rotate(180,%s, %s)" % (x, y) # Set text rotation to upside down
syn_rot = "rotate(0,%s, %s)" % (x, y) # Set text rotation to normal
txt_rotate = "rotate(270,%s, %s)" % (x, y) # Set label rotation to vertical
# Symbolize different line types (anticlines, synclines, etc.)
if str(line_type)[:4].lower() == "anti":
geolines_layer.add(dwg.polyline(line_points, stroke="magenta", stroke_width=0.7087))
# Label the tick "U" (upside down)
geolines_layer.add(dwg.text("U", insert=(x, y + 10),
font_size="6px",
fill="magenta",
text_anchor="middle",
transform=anti_rot))
if line_label:
geolines_layer.add(dwg.text(line_label, insert=(x + 5, y + 5),
font_size="6px",
fill="black",
transform=txt_rotate))
elif str(line_type)[:3].lower() == "syn":
geolines_layer.add(dwg.polyline(line_points, stroke="magenta", stroke_width=0.7087))
# Label the tick "U" (right side up)
geolines_layer.add(dwg.text("U", insert=(x, y - 5),
font_size="6px",
fill="magenta",
text_anchor="middle",
transform=syn_rot))
if line_label:
geolines_layer.add(dwg.text(line_label, insert=(x + 5, y + 5),
font_size="6px",
fill="black",
transform=txt_rotate))
else:
geolines_layer.add(dwg.polyline(line_points, stroke="yellow", stroke_width=0.5))
# Only label the tick with line_type if it is not a fold
if line_label:
geolines_layer.add(dwg.text(line_type + "; " + line_label, insert=(x + 5, y + 1.5),
font_size="6px",
fill="black",
transform=txt_rotate))
else:
geolines_layer.add(dwg.text(line_type, insert=(x + 5, y + 1.5),
font_size="6px",
fill="black",
transform=txt_rotate))
addMsgAndPrint(f" GEL intersections written to SVG file...")
# Process line intersections for MUL tick marks
if len(mulines_fc) > 3:
with arcpy.da.SearchCursor(mulines_fc, ["SHAPE@", "Label"]) as cursor:
for row in cursor:
line_geom, line_label = row
intersection_points = polyline.intersect(line_geom, 1)
for point in intersection_points:
"""addMsgAndPrint(f"Map coords of MUL-profile intersection: X: {point.X}, Y: {point.Y}")"""
dist = calculate_distance(polyline.firstPoint, point)
x, y = map_to_svg_coords(dist, point.Z) # Location on 2D profile
# Define the tick as a polyline with these vertices
line_points = [(x, y + 3), (x, y), (x, y - 3)]
# Symbolize line and label
geolines_layer.add(dwg.polyline(line_points, stroke="lime", stroke_width=0.5669))
txt_rotate = "rotate(0,%s, %s)" % (x, y) # Set text rotation to normal
if line_label:
geolines_layer.add(dwg.text(line_label, insert=(x, y - 5),
font_size="6px",
font_family="FGDCGeoAge",
fill="lime",
text_anchor="middle",
transform=txt_rotate))
addMsgAndPrint(f" MUL intersections and labels written to SVG file...")
# Process map unit polygon labels along the profile
if len(mupolys_fc) > 3:
with arcpy.da.SearchCursor(mupolys_fc, ["SHAPE@", "Label"]) as cursor:
for row in cursor:
polygon_geom, polygon_label = row
intersection_line = polyline.intersect(polygon_geom, 2)
if intersection_line:
# Get the centroid of the line
label_point_int_line = intersection_line.labelPoint
# Check the line for multiple parts which occurs when a polygon overlaps the line more than once
if intersection_line.isMultipart and intersection_line.partCount > 0:
"""addMsgAndPrint(f"intersection_line part count for {polygon_label}: {intersection_line.partCount}")"""
for i in range(intersection_line.partCount):
# Get the array of points that make up the part
this_part = intersection_line.getPart(i)
# Convert the array into a polyline
part_polyline = arcpy.Polyline(this_part)
# Calculate the distance from the beginning of the profile to the beginning of the part
dist_start = calculate_distance(polyline.firstPoint, part_polyline.firstPoint)
# Calculate the distance from the beginning of the profile to the end of the part
dist_end = calculate_distance(polyline.firstPoint, part_polyline.lastPoint)
# Calculate the distance from the beginning of the profile to the middle of the part
dist_middle = dist_start + ((dist_end - dist_start) / 2)
# Determine the X coordinate label position
x_start, _ = map_to_svg_coords(dist_start, min_z)
x_end, _ = map_to_svg_coords(dist_end, min_z)
x_label = (x_start + x_end) / 2
# Add y_label position based on the closest Z coord to the middle of the part
mup_point_z = find_closest_Z_value(profile_points, dist_middle)
_, y_label = map_to_svg_coords(dist_start, mup_point_z)
"""addMsgAndPrint(
f" part {i} for {polygon_label} with {len(this_part)} points plot label at {x_label}, {y_label}")"""
if polygon_label:
mupolys_layer.add(
dwg.text(polygon_label,
insert=(x_label, y_label - 3.5),
font_family="FGDCGeoAge",
font_size="6px",
fill="black",
text_anchor="middle"))
else:
dist_start = calculate_distance(polyline.firstPoint, intersection_line.firstPoint)
dist_end = calculate_distance(polyline.firstPoint, intersection_line.lastPoint)
x_start, _ = map_to_svg_coords(dist_start, min_z)
x_end, _ = map_to_svg_coords(dist_end, min_z)
x_label = (x_start + x_end) / 2
# Add y_label position based on Z coord of label point of intersection_line segment
_, y_label = map_to_svg_coords(dist_start, label_point_int_line.Z)
"""addMsgAndPrint(
f"intersection_line is single part for {polygon_label} with {intersection_line.pointCount} points plot label at {x_label}, {y_label}")"""
if polygon_label:
mupolys_layer.add(
dwg.text(polygon_label,
insert=(x_label, y_label - 3.5),
font_family="FGDCGeoAge",
font_size="6px",
fill="black",
text_anchor="middle"))
addMsgAndPrint(f" MUP labels written to SVG file...")
# Select and plot nearby points as ball-and-stick symbols
if len(orpoints_fc) > 3:
arcpy.MakeFeatureLayer_management(orpoints_fc, "temp_points_layer")
arcpy.SelectLayerByLocation_management("temp_points_layer", "WITHIN_A_DISTANCE", polyline, search_distance)
with arcpy.da.SearchCursor("temp_points_layer", ["SHAPE@", "Type", "Azimuth", "Inclination"]) as cursor:
for row in cursor:
point, fType, azimuth, inclination = row
# Only process points that have the attributes we're looking for
if (inclination != -9999) and (str(fType)[:7] == "bedding"):
# Get the distance along the profile of the projected location of the selected point
dist = sqrt((point.getPart().X - polyline.firstPoint.X) ** 2 + (
point.getPart().Y - polyline.firstPoint.Y) ** 2)
# Finds the Z value of the closest "dist" profile point
orp_point_z = find_closest_Z_value(profile_points, dist)
# Get plot x, y of the ball based on the distance along the line and the Z value of the profile
x_ball, y_ball = map_to_svg_coords(dist, orp_point_z)
# Calculate the location of the end of the stick using trig
if azimuth < 180:
# Add 90 degrees to inclination value to convert to Arithmetic rotation (dip/plunge is east)
y_tip = y_ball - 10 * cos(radians(inclination + 90))
x_tip = x_ball + 10 * sin(radians(inclination + 90))
else:
# Reverse the x direction value to simulate dip/plunge west
y_tip = y_ball - 10 * cos(radians(inclination + 90))
x_tip = x_ball - 10 * sin(radians(inclination + 90))
# Draw the ball-and-stick symbol
dip_layer.add(dwg.circle(center=(x_ball, y_ball), r=1.5, fill="green"))
dip_layer.add(
dwg.line(start=(x_ball, y_ball), end=(x_tip, y_tip), stroke="green", stroke_width=0.5))
addMsgAndPrint(f" ORP data projected to profile line and written to SVG file...")
# Save the SVG file
#dwg.save()
# Save the SVG file
dwg.save()
addMsgAndPrint(f"SVG file saved to {output_svg}.")
### COLLECT USER INPUTS ###
# Get start time:
start_time = datetime.datetime.now()
# Set workspace environment
#arcpy.env.workspace = r"C:\temp\geodatabase" or r"C:\temp\folder"
output_ws = sys.argv[1]
arcpy.env.workspace = output_ws
# Input polyline feature class (gdb fc or shapefile)
#input_fc = "input_polyline_feature_class"
input_fc = sys.argv[2]
# DEM raster; must be same projection as input_fc, z units must match x, y units or use VE to adjust
#dem = r"C:\temp\dem.tif"
dem = sys.argv[3]
# Sample interval in map units
sample_interval = sys.argv[4]
# Output PolylineZM feature class name (gdb fc or shapefile)
#output_fc_name = "output_polylineZM_feature_class"
output_fc_name = sys.argv[5]
# Output SVG file path
#output_svg = r"C:\temp\folder\svg_file.svg"
output_svg = sys.argv[6]
# Interval for length grid (in projection units)
#interval_length = 100
interval_length = sys.argv[7]
# Interval for elevation grid (in elevation units)
#interval_elevation = 10
interval_elevation = sys.argv[8]
# Map scale (scale denominator)
#map_scale = 24000
map_scale = sys.argv[9]
# Vertical Exaggeration (default is 1)
#vertical_ex = 1
vertical_ex = sys.argv[10]
### Additional inputs to plot on profile:
# ticks_and_anno = True
ticks_and_anno = sys.argv[11]
# ContactsAndFaults
#CAF_fc = r"path_to_your_CAF_feature_class"
CAF_fc = sys.argv[12]
# GeologicLines
#CAF_fc = r"path_to_your_GEL_feature_class"
GEL_fc = sys.argv[13]
# MapUnitLines
#MUP_fc = r"path_to_your_MUL_feature_class"
MUL_fc = sys.argv[14]
# MapUnitPolys
#MUP_fc = r"path_to_your_MUP_feature_class"
MUP_fc = sys.argv[15]
# OrientationPoints
#ORP_fc = r"path_to_your_ORP_feature_class"
ORP_fc = sys.argv[16]
# Search (buffer) distance along profile line to look for OrientationPoints
#search_distance = 50 # Search distance for nearby points
search_distance = sys.argv[17]
### RUN SUBS TO CREATE OUTPUTS ###
# Create PolylineZM feature class path from workspace and class name
#feature_class = r"path_to_your_polyline_feature_class"
feature_class = os.path.join(output_ws, output_fc_name).replace("\\", "/")
spatial_ref = arcpy.Describe(input_fc).spatialReference
# Process the polyline and create the PolylineZM
process_polyline_zm(input_fc, dem, feature_class, spatial_ref)
# Analyze ZM polyline feature class and create SVG output
plot_elevation_profile(feature_class, output_svg, interval_length, interval_elevation, map_scale, vertical_ex, CAF_fc,
GEL_fc, MUL_fc, MUP_fc, ORP_fc)
# Inform the user how long the process took
end_time = datetime.datetime.now()
delta = end_time - start_time
addMsgAndPrint(f"Processing time: {delta.total_seconds()} seconds")
reacted with thumbs up emoji reacted with thumbs down emoji reacted with laugh emoji reacted with hooray emoji reacted with confused emoji reacted with heart emoji reacted with rocket emoji reacted with eyes emoji
-
If you or someone in your organization still prefers to draw cross sections "the old fashioned way" then this script is for you!
Supply a feature class with a single polyline and specify a DEM (GRID, tiff, img, etc.; NOT compatible with mosaic datasets... yet) and it generates a PolylineZ (3D) feature class and a 2D profile plot and exports it to an SVG file. Note that the spatial references of all inputs must match and the script assumes that the x, y, and z values are in meters. You can specify the horizontal and vertical grid spacing of the plot as well as the vertical exaggeration and map scale. Optionally annotates the profile with info from the user-specified GeMS feature classes. The SVG file is written using the Inkscape SVG format which allows the use of "layers" that Illustrator interprets as groups, keeping the grid, profile, and annotations organized and easy to separate.
I had created a basic profile graphic tool many years ago to solve the problem that the ArcMap 3D/Spatial Analyst profile tool produced an accurate but unsatisfactory profile graphic. The ArcMap output was not scaled to the map layout and the horizontal and vertical scales were in no way related and none of these options could be specified which required scaling the plot later in Illustrator. Then the geologic map information had to be added manually by superimposing the profile on the map and transferring the data to the profile. My script created a scaled profile, but it did not add any geologic information as the lack of a common schema to the inputs made this difficult to automate.
After a few fits and starts, I finally finished this script when I had to produce new profiles for one of our geologists to complete some cross sections for a previously unpublished map that I was compiling. Drafts of the cross sections had been hand drawn, but the profiles were generalized a bit too much for our liking and the contacts were adjusted during the compilation process making the hand drawn cross sections obsolete.
Key features of this script:
Check out the code below for a detailed description and see an example of the output in the attached SVG file. The example is a re-creation of the profile shown in cross section B-B' on the Bacon Gap geologic quadrangle map: https://ngmdb.usgs.gov/Prodesc/proddesc_108498.htm. The cross section on the published map was created by the author; I only re-created the profile and plotted the pertinent geologic map data to show how the tool output can be used as the starting point for drawing a cross section.
I have also attached an ArcGIS Pro .atbx file with the script in a zip file to make setting up the tool easier since it has a lot of parameters. Unzip, add the .atbx file to your ArcGIS Pro Favorites, verify in the script properties that the path to the .py file is correct, and run the tool. The tool does have one additional Python dependency (svgwrite) that you may have to add to your Python environment if you don't already have it. In the folder of your cloned ArcGIS Python environment, install svgwrite (Start > Python command prompt > pip install svgwrite) and the code should work when run from an ArcGIS session.
GeneratePolylineZExport2DProfileSVG.zip
Beta Was this translation helpful? Give feedback.
All reactions