Topo to Raster Batch Process question (Python newbie)

[AndrewMaracini]

Hi,

I've searched and researched and sampled but I'm hitting the wall here.

I'm trying to run Topo to Raster on 500+ inputs (topo contours and lake boundaries). I've tried to follow samples and examples in books and I think my workflow is going in the right direction but I'm getting an error, and hoping that someone can find my mistake. I've got all of the files in one directory and I defined my variables with a wildcard which is where I think my problem is but I'm not sure how to define the variable. For example the topo line feature is made up of  519 seperate files that have the naming convention "topo_ras_input_topoSection1","topo_ras_input_topoSection2" The lakes boundary file have a similar naming convention: "topo_ras_input_watersecSection1".

After I defined my variables I then list all of the feature classes and then use a for statement to start iterating through the list, I really should end up with 519 rasters as output but I'm not sure how to specify that either.

My code (I use that term loosely) is attached. Thanks in advance for your help.

[MarcinGasior]

Assuming that all your shapefiles are in one catalog and the follow naming convention "topo_ras_input_topoSection1","topo_ras_input_topoSection2", ... , "topo_ras_input_topoSection519", the following script should work:

import arcpy
arcpy.env.workspace = "D:/GIS_data/Topo_Water"
arcpy.CheckOutExtension("Spatial")

cellSize = 30 #adjust this value to your data

#assume there's 519 sections in workspace(you can test script with smaller value, eg. 3)
for i in range(1, 520):
    myTopoContour = arcpy.sa.TopoContour([['topo_ras_input_topoSection'+str(i) + '.shp','contour']])
    myTopoLake = arcpy.sa.TopoLake(['topo_ras_input_watersecSection'+str(i) + '.shp'])

    outTopoToRaster = arcpy.sa.TopoToRaster([myTopoContour, myTopoLake], cellSize)
    outTopoToRaster.save("D:/GIS_data/Raster/twraster" + str(i).zfill(3))

The output raster name follow the rule: twraster001, twraster002, ...
You can try running the script for testing for few shapefiles - adjust number in range(1, x) part.
You may need to test to find optimal cell size of output raster.

[markdenil]

Try being explict with setting the names.
Even if inputing the bald wildcards worked, you would get one output from all 519 input pairs, not 519 outputs...
This is not tested, but may point you in the right direction.
It assumes there are fc and raster pairs numbered the same.
It itterates through the feature classes, assuming there is a matching raster.
It uses the same 'contour' fc input for each loop.
It outputs a twraster<number> where <number> matches the input number.

# ---------------------------------------------------------------------------
# Listdata.py
# # Description: 
# Lists data sets in a workspace and then runs Topo to Raster on each feature
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
# Check out any necessary licenses
arcpy.CheckOutExtension ("Spatial")
# Set the workspace for the ListFeatureClass function
arcpy.env.workspace = "D:/GIS_data/Topo_Water"
# Use the ListFeatureClasses function to return a list of 
#  all shapefiles and feature classes.
fcList = arcpy.ListFeatureClasses()

# Process each file with loop in topo to raster
for fc in fcList:
    fcBase = fc.rstrip(".shp")  #strips off the '.shp.'
    fcNUM = fcBase .lstrip("topo_ras_input_topoSection")  # should leave just the digits
    rasIN = 'topo_ras_input_watersecSection%s' % (str(fcNUM)) # output for this iteration
 
    myTopoContour = TopoContour([[fc,'contour']])
    myTopoLake = TopoLake([rasIN]) # does not check if this exists....

    # Execute TopoToRaster
    outTopoToRaster = TopoToRaster([myTopoContour, myTopoLake])
# Save the output 
    outTopoToRaster.save("D:/GIS_data/Raster/twraster%s" % (fcNUM)) # saves each one

... and I see that as I was typing, someone else answered too....

[AndrewMaracini]

Assuming that all your shapefiles are in one catalog and the follow naming convention "topo_ras_input_topoSection1","topo_ras_input_topoSection2", ... , "topo_ras_input_topoSection519", the following script should work:

import arcpy
arcpy.env.workspace = "D:/GIS_data/Topo_Water"
arcpy.CheckOutExtension("Spatial")

cellSize = 30 #adjust this value to your data

#assume there's 519 sections in workspace(you can test script with smaller value, eg. 3)
for i in range(1, 520):
    myTopoContour = arcpy.sa.TopoContour([['topo_ras_input_topoSection'+str(i) + '.shp','contour']])
    myTopoLake = arcpy.sa.TopoLake(['topo_ras_input_watersecSection'+str(i) + '.shp'])

    outTopoToRaster = arcpy.sa.TopoToRaster([myTopoContour, myTopoLake], cellSize)
    outTopoToRaster.save("D:/GIS_data/Raster/twraster" + str(i).zfill(3))

The output raster name follow the rule: twraster001, twraster002, ...
You can try running the script for testing for few shapefiles - adjust number in range(1, x) part.
You may need to test to find optimal cell size of output raster.

Marcin,

Thanks for your quick reply. Initially the code didn't work but I found that I had an error in my filename. The code runs, thank you, thank you! Can you please explain why there was no need to list features and loop? I thought that was the proper way to approach the problem? Is your approach something that can be used in other cases where I have many files in the same directory?

[MarcinGasior]

Can you please explain why there was no need to list features and loop? I thought that was the proper way to approach the problem? Is your approach something that can be used in other cases where I have many files in the same directory?

Usually creating a list of feature classes with ListFeatureClasses() and looping through this list is good choise. As you notice, Mark's solution uses this approach.
But this is most convenient when you do something with one feature class at once.
In your case you need to use two feature classess and the naming convention allowed to use numeric iterator to loop through the data.

[AndrewMaracini]

Mark,


I've got a follow-up question. I was able to get Marcin's code to run but in doing so found an error in my logic and I think that I need to include an if elif statement in the code in order for the TopoToRaster function to work properly. And I think your code is the better starting point. Because of my lack of experience in programming, I'm in over my head here. I'll try to explain in detail and include a screen capture for good measure.

I'm trying to run TopotoRaster on an entire county at high cell resolution (2ft). In order to do so I had to break the 2ft countours out and water boundary polygon layer into 1 square mile sections. Since not all sections have water and not all sections have topo I had to remove all shapefiles with no topo. What I need the program to do is to compare file names to check if the section numbers are equal or even exist and run Topo to Raster with the appropriate parameters such that:

If topoSection == watersecSection then run TopoToRaster([myTopoContour, myTopoLake], cellSize 2)

elif run TopoToRaster([myTopoContour], cellSize2)

I basically need to evaluate the numbers in each filename and if they are equal to include the TopoLake topo class and if not just run the contour. In my screenshot you will see topoSection119 and watersecSection119 in this case run both parameters. In other cases I have topoSection0 and no watersecSection0 where I just need to run one parameter.

I really have no clue how to write the code. I'm running the TopoToRaster Iterator in model builder without the Lake boundary function just to see how it turns out using only the Topo line data, and it seems to be working. If I could use the iterator with the two files I would have done that but I can't figure out how to make that work either.

Thanks in advance,

Andy

Try being explict with setting the names.
Even if inputing the bald wildcards worked, you would get one output from all 519 input pairs, not 519 outputs...
This is not tested, but may point you in the right direction.
It assumes there are fc and raster pairs numbered the same.
It itterates through the feature classes, assuming there is a matching raster.
It uses the same 'contour' fc input for each loop.
It outputs a twraster<number> where <number> matches the input number.

# ---------------------------------------------------------------------------
# Listdata.py
# # Description: 
# Lists data sets in a workspace and then runs Topo to Raster on each feature
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
# Check out any necessary licenses
arcpy.CheckOutExtension ("Spatial")
# Set the workspace for the ListFeatureClass function
arcpy.env.workspace = "D:/GIS_data/Topo_Water"
# Use the ListFeatureClasses function to return a list of 
#  all shapefiles and feature classes.
fcList = arcpy.ListFeatureClasses()

# Process each file with loop in topo to raster
for fc in fcList:
    fcBase = fc.rstrip(".shp")  #strips off the '.shp.'
    fcNUM = fcBase .lstrip("topo_ras_input_topoSection")  # should leave just the digits
    rasIN = 'topo_ras_input_watersecSection%s' % (str(fcNUM)) # output for this iteration
 
    myTopoContour = TopoContour([[fc,'contour']])
    myTopoLake = TopoLake([rasIN]) # does not check if this exists....

    # Execute TopoToRaster
    outTopoToRaster = TopoToRaster([myTopoContour, myTopoLake])
# Save the output 
    outTopoToRaster.save("D:/GIS_data/Raster/twraster%s" % (fcNUM)) # saves each one

... and I see that as I was typing, someone else answered too....

[markdenil]

It sounds like you have two 'optional if present' inputs to the process.
You want to run the process on ALL the topos present, although some of the tiles may be missing (so there would be nothing to run)
If there is a water layer present, you would want to include it.

I recall that the topos are features and the waters are rasters

The code I suggested already cycles through the vector features, so it naturally weeds out missing tiles.
I would test for the existance of the water raster immidietly after the name is set

myTopoLake= 'topo_ras_input_watersecSection%s' % (str(fcNUM)) 
waterBool = arcpy.Exists(myTopoLake)

Then use the existance boolean to decide which version of the command you need, with or without water.

if waterBool:
    arcpy.TopoToRaster_sa([myTopoContour, myTopoLake], 2)
else:
    arcpy.TopoToRaster_sa([myTopoContour], 2)

[AndrewMaracini]

Mark,

The water layer is also a vector (polygon) shapefile, do I need to modify that code for the water layer? Thanks, for the help!

andy

It sounds like you have two 'optional if present' inputs to the process.
You want to run the process on ALL the topos present, although some of the tiles may be missing (so there would be nothing to run)
If there is a water layer present, you would want to include it.

I recall that the topos are features and the waters are rasters

The code I suggested already cycles through the vector features, so it naturally weeds out missing tiles.
I would test for the existance of the water raster immidietly after the name is set

myTopoLake= 'topo_ras_input_watersecSection%s' % (str(fcNUM)) 
waterBool = arcpy.Exists(myTopoLake)

Then use the existance boolean to decide which version of the command you need, with or without water.

if waterBool:
    arcpy.TopoToRaster_sa([myTopoContour, myTopoLake], 2)
else:
    arcpy.TopoToRaster_sa([myTopoContour], 2)

[markdenil]

The only issue is that you are looping through a list of all the feature classes.
If that list includes both the topo and the water, then things might get confused.

Use a wildcard on ListFeatureClasses so you only get the contours in the list.
fcList = arcpy.ListFeatureClasses("topo_ras_input_topoSection*")
Otherwise, there shouldn't be an issue with finding the right pairs of inputs.

I would try a version of the process that does not do the actual TopoToRaster,
but just prints out the names of the inputs instead.
You can see if you are getting the right files for each itteration.
something like:

if waterBool:
    print "%s     %s" % (myTopoContour, myTopoLake)
else:
     print "%s" % (myTopoContour)

[AndrewMaracini]

Thanks Mark,

I'm getting a NameError that I can't seem to fix. I've tried changing indents, renaming but still keep getting the following error:

Traceback (most recent call last):
File "D:/Python_code/andy_scripts/DenilCode_AJM_test.py", line 24, in <module>
myTopoContour = TopoContour([[fc,'contour']])
NameError: name 'TopoContour' is not defined

My code is attached. Note that the filenames have changed since the original code, I dropped the original prefix.

The only issue is that you are looping through a list of all the feature classes. 
If that list includes both the topo and the water, then things might get confused. 

Use a wildcard on ListFeatureClasses so you only get the contours in the list. 
fcList = arcpy.ListFeatureClasses("topo_ras_input_topoSection*") 
Otherwise, there shouldn't be an issue with finding the right pairs of inputs. 

I would try a version of the process that does not do the actual TopoToRaster,  
but just prints out the names of the inputs instead. 
You can see if you are getting the right files for each itteration. 
something like: 

if waterBool:
    print "%s     %s" % (myTopoContour, myTopoLake)
else:
     print "%s" % (myTopoContour)