Saturday 23 May 2015

spatial analyst - Using Loop with Raster Calculator in ArcPy?


I am trying to loop through multiple rasters and extract cells that are above 0.1 The code I am trying to use is:


import arcpy 
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'F:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master'
#pathway to all rasters in workspace directory

rasters = arcpy.ListRasters('*.tif*')
for raster in rasters:
#give new file a name
outraster = raster.replace('.tif','_water.tif')
arcpy.gp.RasterCalculator_sa("""raster >= 0.1""",outraster)
print('Done Processing')

For some reason I can't copy and paste the error (python actually shuts down when I run this code) but here is a screen shot of part of it.


enter image description here


EDIT: My parameters have changed and I now need everything >=-0.4 and but !=0.00. I am trying to use:



import arcpy 
from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'D:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master'
#pathway to all rasters in workspace directory
rasters = arcpy.ListRasters('*.tif*')
for raster in rasters:
#Save temp raster to disk with new name
ras = Raster(raster)
outraster = Con(ras >= -0.4 & ras !=0.00, ras)

outraster.save(raster.replace('.tif','_water.tif'))
print('Done Processing')

but this returns:


ValueError: The truth value of a raster is ambiguous. Invalid use of raster with Boolean operator or function. Check the use of parentheses where applicable.

EDIT:


I think I just needed to change this line:


 outraster = Con((ras >= -0.4) & (ras !=0.00), ras)

Answer




From the ArcGIS Help:



The Raster Calculator tool is intended for use in the ArcGIS Desktop application only as a GP tool dialog box or in ModelBuilder. It is not intended for use in scripting and is not available in the ArcPy Spatial Analyst module.



You should use arcpy.sa map algebra instead:


import arcpy 
from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'F:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master'
#pathway to all rasters in workspace directory

rasters = arcpy.ListRasters('*.tif*')
for raster in rasters:
outraster = Raster(raster) >= 0.1
#Save temp raster to disk with new name
outraster.save(raster.replace('.tif','_water.tif'))
print('Done Processing')

Note: both your original expression and the map algebra version above will not extract the values of the cells, but output a 1 where the condition is met and 0 otherwise. If you want the actual cell values, you need to use the Con function:


ras = Raster(raster)
outraster = Con(ras >= 0.1, ras)


You may be having issues with the spaces in your path, it would be better to remove the spaces, but you could also try setting the scratch workspace to a folder without spaces, i.e arcpy.env.scratchWorkspace = r'C:\Temp'


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...