Python script coordinates system problem

Discussion created by teabox on Jun 24, 2011
Latest reply on Jun 28, 2011 by l3ll0
Hi !

I wrote the following python script and I ran it one time and all works perfectly. Then I changed the name of the output folder in order to create an other surface for and other discharge. At the end all the output files are correct excepted the output file of the SetNull function. This file is created with no coordinates system and the value of the raster are "consequently" wrong. I tried several time to rewrite the script but I still have the same problem. I do not understand why the coordinates system is not set up in this only file when all the other output files are correct.
Could someone help me ?



# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
env.workspace = "C:/Users/Teabox/Documents/ArcGIS/MasterThesis" #Set your environment file
arcpy.env.outputCoordinateSystem = "D:/Programmes/ArcGisDesktop10/Desktop10.0/Coordinate Systems/Projected Coordinate Systems/UTM/WGS 1984/Northern Hemisphere/WGS 1984 UTM Zone 32N.prj"
arcpy.env.cellSize = "C:/Users/Teabox/Documents/ArcGIS/MasterThesis/lundes_raster"

# Set local variables
inPntFeat = "EndPoints.shp" #Set the name of the input file
zField = "elevation" #Set the name of the interpolate field
cellSize = arcpy.env.cellSize = "C:/Users/Teabox/Documents/ArcGIS/MasterThesis/lundes_raster"
splineType = "REGULARIZED"
weight = 0.1

# Check out the ArcGIS Spatial Analyst extension license

print "Execute Spline"
# Execute Spline
outSpline = Spline(inPntFeat, zField, cellSize, splineType, weight)

# Save the output

print "Execute Raster Calculator"
#Evaluate the difference between the interpolated water surface and the digital elevation model of the river bed
#The difference gives the areas covered by water as positive number and the rest as negative numbers.
outRas = Raster("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/LowWaterL") - Raster("lundes_raster")

# Final step to display only the values above zero to create a final layer showing only the area covered by water.
outSetNull = SetNull("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/output_raster", "C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/output_raster", "VALUE < 0")

print "Reclassify water level layer"
# Reclassify the water level layer in order to run the ZonalGeometry tools.
# This Zonal tools works only on integer raster.
Waterlevel = "C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/area.img"
outReclassWL = Reclassify( Waterlevel, "Value", RemapRange([[0,1,1], [1,3,1]]), "NODATA")  

print "Run the ZonalGeometry tool"

#ZonalGeometry tool
# Set local variables
inZoneData = ("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/ReclassWL")
zoneField = "Value"

# Execute ZonalStatistics
outZonalGeometry = ZonalGeometry(inZoneData, zoneField, "AREA", cellSize) 

# Save the output

#Extract the value of the water area
WaterArea = arcpy.GetRasterProperties_management("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/zonegeo", "MAXIMUM")
print "Total water surface in m²"
print WaterArea

print "End"