1

I have created profiles at multiple locations with lines.

Now I want to add points to that graph, how do I do that?

Here is the script:

import arcpy
import os
import matplotlib.pyplot as plt
import numpy

#settings for datasources
linesPath = "lines"
rasterPath = "raster"
pointsPath = "points"

# settings for output profiles PNG's
profileFolder = r"\test_arcpy\profile"
profilePrefix = "ArcPy_"
profilePostfix = "_new.png"
fldName = "profile_na"

# standard NoData mapping
NoDataValue = -9999

# describe raster
inDesc = arcpy.Describe(rasterPath)
rasMeanCellHeight = inDesc.MeanCellHeight
rasMeanCellWidth = inDesc.MeanCellWidth

# search cursor
fldShape = "SHAPE@"
flds = [fldShape, fldName]
with arcpy.da.SearchCursor(linesPath, flds) as curs:
    for row in curs:
        feat = row[0]
        profileName = row[1]

        # read extent of feature and create point for lower left corner
        extent = feat.extent
        print "Processing: {0}".format(profileName)
        pntLL = arcpy.Point(extent.XMin, extent.YMin)

        # determine number of rows and cols to extract from raster for this feature
        width = extent.width
        height = extent.height
        cols = int(width / rasMeanCellWidth)
        rows = int(height / rasMeanCellHeight)

        # create Numpy array for extent of feature
        arrNP = arcpy.RasterToNumPyArray(rasterPath, pntLL, cols, rows, NoDataValue)

        # create empty arrays for distance (d) and altitude (z) and NAP line (zv)
        d = []
        z = []
        zv = []

        # loop through polyline and extract a point each meter
        for dist in range(0, int(feat.getLength("PLANAR"))):
            XYpointGeom = feat.positionAlongLine(dist, False)
            XYpoint = XYpointGeom.firstPoint

            # translate XYpoint to  row, col (needs additional checking)
            c = int((XYpoint.X - pntLL.X) / rasMeanCellWidth)
            r2 = int((XYpoint.Y - pntLL.Y) / rasMeanCellHeight)
            r = rows - r2
            if c >= cols:
                c = cols-1
            if r >= rows:
                r = rows-1

            # extract value from raster and handle NoData
            zVal = arrNP[r,c]
            if not zVal == NoDataValue:
                d.append(dist)
                z.append(arrNP[r,c])
                zv.append(0)

        # define output profile filename and figure size settings
        elevPNG = profileFolder + os.sep + profilePrefix + profileName + profilePostfix
        fig = plt.figure(figsize=(10, 3.5)) # inches

        # plot profile, define styles
        plt.plot(d,z,'r',linewidth=0.75)
        plt.plot(d,z,'ro',alpha=0.3, markersize=3)
        plt.plot(d,zv,'k--',linewidth=0.5)
        plt.xlabel('Distance from start')
        plt.ylabel('Elevation')
        plt.title('Profile {0} using Python matplotlib'.format(profileName))

        # change font size
        plt.rcParams.update({'font.size': 8})

        # save graph as PNG
        fig.savefig(elevPNG, dpi=300)

# delete some objects
del arrNP, XYpointGeom, XYpoint
print "PNGs stored in: {0}".format(profileFolder)

enter image description here

  • At the moment it looks like you are presenting a lot of code, with little description, and asking us to improve it. I suspect that is why no potential answerers have volunteered their time to try and answer it. – PolyGeo Aug 2 '19 at 13:07

Browse other questions tagged or ask your own question.