Need python equivalent to Avenue’s “ReturnContour”

[RebeccaStrauch__GISP]

I’m working on updating one of my old Avenue scripts.

Short description (of this particular problem in a longer process): given a set of random points within a study area, I need to create a contour at the elevation value of the random point, and using that as the center, create a transect on that contour at a given length, 50% in each direction….adjusting if necessary to keep within the study area.  (and there are other factors I don’t need to get into right now).

I’ve been able to grab the xy from a given point, and retrieve the value from the raster. The problem I am having right now is I can’t find an arcpy replacement for  ReturnContour (as shown in this line)

theFullContour = ElevTheme.GetGrid.ReturnContour(p, z, Prj.MakeNull)

My guess is it a NumPy solution (Dan Patterson maybe??), but I’m not sure what to search for.  Once I get the (poly)line, I should be ok with geometries.)

Does anyone have a simple way in python to return a single contour given an x/y/z and a raster?  Anyone in Geoprocessing or Imagery and Rasters  ??  Please, python solution only....no C#/arcobjects.

Thanks!

EDIT:  meant to include a snippet of test code for getting the point, in case anyone wants to test on their own data.  I'm grabbing point 7 of 2000 because I know it will return a good contour (some will not....too flat, which is handled elsewhere)

myPoints = r"C:\ADFG.gdb\myRandomPts2"
theDEM = r"C:\ADFG.gdb\DEM20aExtent"
cnt = 1
for row in arcpy.da.SearchCursor(myPoints, ["SHAPE@XY"]):
    if cnt == 7:
        x,y = row[0]      
        print("{}, {}, {}".format(cnt, x, y))
    else:
        pass
    cnt = cnt + 1
    
#x, y = 377906.9379, 1538997.6698
print("{}, {}".format(x, y))
ptLoc = str(x) + ' ' + str(y)

z1 = arcpy.GetCellValue_management("DEM20aExtent",ptLoc)
print(z1)

Message was edited by: Rebecca Strauch added code snippet

[DanPatterson_Retired · Accepted solution]

Rebecca... Tooootally not checked out, but on the list.  One of the problems is that the Spatial Analyst tools and arcpy.sa don't provide much in the way of interacting with the dataframe, however, this...

Contour List—Help | ArcGIS for Desktop holds some promise

ContourList (in_raster, out_polyline_features, contour_values)

Parameter	Explanation	Data Type
in_raster	The input surface raster.	Raster Layer
out_polyline_features	The output contour polyline features.	Feature Class
contour_values [contour_value,...]	List of z-values for which to create contours.	Double

​IF you can pass the z-value of the location to the contour_values list.  Totally conjecture, but it seems to be the only on that doesn't require you to contour the whole surface specifying a base contour and interval.  Give it a shot

EDIT

just saw your edit, see if you can throw that z value into the contour list

I would suggest being specific rather than general  when using the tool

arcpy.sa.ContourList(... , ... , [elevation value in list] )

[RebeccaStrauch__GISP]

Thanks Dan.  I had noticed that command, but didn't look at it close enough, but it looks promising.  I added

outCountor = r"C:\ADFG.gdb\c1"
from arcpy.sa import *
ContourList(theDEM, outCountor, z1)

And it hasn't crashed....might be an overnight run.  But if that works, I'll have to work on getting the extent to be much smaller that it is right now (about 12,100 sq mi...110x110)

Thanks for the fast feedback....this really is a great community.

[DanPatterson_Retired]

I hope z1 was in a list... but sometimes esri tools require a list unless only one value is given  (PS always test on a small sample )

[RebeccaStrauch__GISP]

z1 was just one point, which was what I need.  Anyway, it just finished.....about 10 minutes to run the one point/elevation.  I'll have to work on speeding it up since I could easily have 2000 or more to run.

It of course returns too much right now, but that again is something that I can work on. 

Thanks Dan, that was exactly the push I needed in the right direction!

[SteveLynch]

Rebecca

Remember to use SnapRaster (=raster being contoured) as well when setting the extent as ContourList will resample the raster before contouring if the extent does not fall exactly on the origin of a cell.

-Steve

[RebeccaStrauch__GISP]

Thanks Steve.  I will give that a try, but that may require me to move the random seed point (x/y) since that is used, (more-or-less, may slide to stay in study area), as the center point of the transect. I'll have to make sure that doesn't negatively impact the intent of the randomness.

I will be processing one random point at a time, and will only need about 20k of the contour (10k either side of point, along contour....adjusted if it hits the study boundary).  I'm thinking about trying to maybe use an in_memory version of the raster, maybe clipping the raster to an extent 40k on either size of the point.  I'm not sure if this is possible, but I would think a smaller raster extent would speed up the process.

Just thinking "out loud"..One thing the ContourList command may help with is locating "continuation" contours, for example, when the contour is < 20k (e.g., close to top of mountain).  In that case, the plane goes off transect, and restarts on another mountain at same elevation. If I can choose a second contour line for them to continue, it can take some of the human influence out of the randomness.  (if that makes sense)    I have most of my logic (except the "continuation" option) already written in my Avenue program that I'm translating.  But I have to say, ReturnContour in AV3.x works much faster than ContourList at least without setting up environments....so I will try your suggestion.

EDIT: Steve, as you mentioned, the SnapRaster really sped up the process (to < 15 seconds), and the contour results are the same.  The point isn't on the contour itself, in either case, but I can adjust that as needed.  On to the next step.  Thanks for that hint!

[SteveLynch]

There are a few other ways to do this 🙂

1. Load the Spatial Analyst toolbar (Customize->Toolbars -> Spatial Analyst). The Create Contour tool will do what you want albeit at screen resolution, i.e. the more you zoom the more detail you'll get, or
2. Run some geoprocessing tools
   1. ExtractValuesToPoints (using the single point and the raster as input) and remember to check "Interpolate values at the point locations"
   2. in the output copy the value from RASTERVALU
   3. ContourList, use the raster and paste the copied value into the Contour values control.
   4. SelectByLocation use the contour output and the single point with relationship=WITHIN_A_DISTANCE and maybe a small search distance

-Steve

BTW, from the Avenue doc

The following example creates a Graphic line of a contour passing through a point specified with the cursor on the active View. The contour is created from an active GTheme. The value of the contour is displayed in the status bar. The example must be executed from an Apply event since it uses ReturnUserPoint to get the point. It also assumes that there is only one Theme active and that it is a GTheme.

theView = av.GetActiveDoc

theGridTheme = theView.GetActiveThemes.Get(0)

thePoint = theView.GetDisplay.ReturnUserPoint

zValue = 0

aPrj = Prj.MakeNull

theLine = theGridTheme.GetGrid.ReturnContour(thePoint, zValue, aPrj)

theGraphicLine = GraphicShape.Make(theLine)

theView.GetGraphics.Add(theGraphicLine)

av.ShowMsg("Contour Value : "++zValue.AsString)

[RebeccaStrauch__GISP]

Thanks Steve,

Many good tips there that I wil check out. 

I actually will possible create (now) five or six different types of transects that depend on the terrain quite a bit (i.e. how super cubs can fly)....one of the new options was first talked about this morning  In the Avenue program I have only three (contour, straight and straight-hinged), but there were a few areas that seemed to get under sampled.  I'm still working on the logic (both from the needs of the biostatistician and program side) since breaking the logic of each step makes life much easier once I start putting all the pieces together.

This has always been a great project to play with new challenges. One of the nice things is I've been here from the start of design and running of this project (1992-ish) so I've very in tune with it. But I wouldn't doubt throw I'll throw some additional questions out here when I run into roadblocks.  It's usually only run every few years, so super efficiency isn't required, but I can see the transect part being used more often by other project if I do it right.

Eventually, I'll have all the script in toolbox and a python addin (not to be confused with a python toolbox). FWIW.