Python script coordinates system problem

[ThibaultBoissy]

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 ?

Regards,

Thibault

# 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
arcpy.CheckOutExtension("Spatial")

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

# Save the output
outSpline.save("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/LowWaterL")

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")
outRas.save("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/output_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")
outSetNull.save("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/area.img")

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")  
outReclassWL.save("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/ReclassWL")

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
outZonalGeometry.save("C:/Users/Teabox/Documents/ArcGIS/MasterThesis/test/zonegeo")

#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"

[curtvprice]

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.

This can happen with map algebra output with some tools, especially mathematical operators like SetNull.

My advice is to copy the coordinate system from your input to the output raster using the Define Projection tool.

[AlessandroCinnirella]

try using


arcpy.env.outputCoordinateSystem = "PROJCS["WGS_1984_UTM_Zone_32N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",9.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",32632]]"

it worked for me with a similar issue.

hope this helps,
AC