Editing Example Python Scripts for Beginners

PaulAlexander1 · ‎02-01-2013

Hi,

I'm a student, I haven't used Python before but am currently working through getting myself up to date on it, I came across a sample python script on the Arc10 Help website that I would like to use.

My understanding is the code batch runs IDW tool on point data that is time enabled, I've tried playing around with the code for 4 days now and am getting no where with it.

Really what I would like to know is what do I have to modify to get it to work with my own data?

I've highlighted the parts of the code (below) that I think I need to modify for it to work with my data but can some one please let me know if there is other parts of the code I need to modify?

Please and thanks in advance:

import arcpy, datetime
# Check out the ArcGIS Spatial Analyst extension for using the IDW interpolation tool
arcpy.CheckOutExtension("spatial")

# Get the layer time properties
lyr = arcpy.mapping.Layer(r"C:\Project\Data\Time\TimeLayer.lyr")
lyrTime = lyr.time

# Calculate the number of iterations based on the time extent and timestep interval
startTime = lyrTime.startTime
endTime = lyrTime.endTime
timeExtent = endTime - startTime
timeStepInterval = lyrTime.timeStepInterval

iterations = timeExtent.days / timeStepInterval.interval

# Get the time field containing the time values associated
# with the data in the time-enabled layer
startTimeField = str(lyrTime.startTimeField)

# Specify the output mosaic dataset to which the interpolated rasters will be added
outputMosaicDataset = r"C:\Project\Output\Output.gdb\outputMosaicDataset"

i = 0
while i <= iterations:
# Formulate the time query and increment the time by the timeStepInterval
currentTime = str(startTime + (i*timeStepInterval))
timeQuery = "\"" + startTimeField + "\"" + " = date '" + currentTime + "'"

# Create an in-memory feature layer containing points that are valid at each timestep
tempFeatureLyr = "tempTimeLayer" + str(i)
arcpy.MakeFeatureLayer_management(lyr, tempFeatureLyr, timeQuery)

# Create an interpolated raster surface using the points valid at each timestep
outRaster = r"C:\Project\Output\Output.gdb\raster" + str(i)
print outRaster
arcpy.gp.Idw_sa(tempFeatureLyr, "Temperature", outRaster)

# Add the newly created raster surface to a Mosaic Dataset
arcpy.AddRastersToMosaicDataset_management(outputMosaicDataset, "Raster Dataset", outRaster)

i = i + 1

# Calculate the statistics on the output Mosaic Dataset for
# classifying your data after new rasters are added
arcpy.CalculateStatistics_management(outputMosaicDataset,"1","1","#")

Zeke · ‎02-01-2013

I don't know the specifics of what you're trying to do, don't work much with rasters, but a couple possible things to look at:

1. Do you have a Spatial Analyst license? Does it actually get checked out? Some error checking would help here. Generically:

try:
   <some code here>
except e as Exception:
   arcpy.AddMessage('some error message ' + e.message

2. 'arcpy.gp.Idw_sa(tempFeatureLyr, "Temperature", outRaster)' doesn't look correct. gp is usually seen in versions of Arc prior to 10, when arcgisscripting was used instead of arcpy. This could be correct, but I've never seen it.
3. It would be helpful if you enclosed your script in code tags. Python is strict about indentation. If there's an error with indentation, we can't see it when the code is just pasted into the post. See the 1st post in this forum on how to paste code.

PaulAlexander1 · ‎02-04-2013

Thanks for your reply,

Here is the code again (below). I should mention again this is from the ArcGIS help site so I don't anticipate any issues with syntax, my novice impression is I simply need to edit the parameters of this script for it to run? Edit: yes I have full license to Arc 10.0 with all extensions.

Specifically what I need to do: I have point data, each point represents a meteorological station and each contains meteorological data which is time enabled (e.g. temperatures for each point at 1am, 2am, 3am etc.) and I was looking for a way to run IDW across all the points on based on the time (so it takes 1am rainfall, runs IDW, then 2am rainfall, runs IDW, then 3am rainfall, run IDW etc.)

Since I'm working with lots of points and more importantly lots of time stamps, I believe the script below is suitable as my understanding is it would sort through these data automatically and run IDW for each time step and then add all the results to a mosaic raster dataset, rather than doing it manually.

I've looked at it again, It seems I need to edit the input layer file to navigate to my own, edit the field that contains my time stamps, edit the output location for the mosaic data set.

I can run the tool and let you know what errors I'm getting back if that would help?

Many thanks again,

Paul.

import arcpy, datetime
# Check out the ArcGIS Spatial Analyst extension for using the IDW interpolation tool
arcpy.CheckOutExtension("spatial")

# Get the layer time properties
lyr = arcpy.mapping.Layer(r"C:\Project\Data\Time\TimeLayer.lyr")
lyrTime = lyr.time

# Calculate the number of iterations based on the time extent and timestep interval
startTime = lyrTime.startTime
endTime = lyrTime.endTime
timeExtent = endTime - startTime
timeStepInterval = lyrTime.timeStepInterval

iterations = timeExtent.days / timeStepInterval.interval

# Get the time field containing the time values associated
# with the data in the time-enabled layer
startTimeField = str(lyrTime.startTimeField)

# Specify the output mosaic dataset to which the interpolated rasters will be added 
outputMosaicDataset = r"C:\Project\Output\Output.gdb\outputMosaicDataset"

i = 0
while i <= iterations:
    # Formulate the time query and increment the time by the timeStepInterval
    currentTime = str(startTime + (i*timeStepInterval))
    timeQuery = "\"" + startTimeField + "\"" + " = date '" + currentTime + "'"

    # Create an in-memory feature layer containing points that are valid at each timestep
    tempFeatureLyr = "tempTimeLayer" + str(i)
    arcpy.MakeFeatureLayer_management(lyr, tempFeatureLyr, timeQuery)

    # Create an interpolated raster surface using the points valid at each timestep
    outRaster = r"C:\Project\Output\Output.gdb\raster" + str(i)
    print outRaster
    arcpy.gp.Idw_sa(tempFeatureLyr, "Temperature", outRaster)

    # Add the newly created raster surface to a Mosaic Dataset
    arcpy.AddRastersToMosaicDataset_management(outputMosaicDataset, "Raster Dataset", outRaster)

    i = i + 1
    
# Calculate the statistics on the output Mosaic Dataset for
# classifying your data after new rasters are added   
arcpy.CalculateStatistics_management(outputMosaicDataset,"1","1","#")

MathewCoyle · ‎02-04-2013

I agree with Greg. I really don't like how the Idw tool is accessed. This should be the proper method.

arcpy.sa.Idw()

It would be helpful if you supplied any error messages you received.

PaulAlexander1 · ‎02-04-2013

Thanks again for your replies,

Ok so here is the message I got back:

Runtime error <type 'exceptions.AttributeError'>: 'Layer' object has no attribute 'time'

after running this code:

import arcpy, datetime
# Check out the ArcGIS Spatial Analyst extension for using the IDW interpolation tool
arcpy.CheckOutExtension("spatial")

# Get the layer time properties
lyr = arcpy.mapping.Layer(r"C:\Users\Paul\Documents\ArcGIS\Python Time\Time Data\Time_Enabled.lyr")
lyrTime = lyr.time

# Calculate the number of iterations based on the time extent and timestep interval
startTime = lyrTime.startTime
endTime = lyrTime.endTime
timeExtent = endTime - startTime
timeStepInterval = lyrTime.timeStepInterval

iterations = timeExtent.days / timeStepInterval.interval

# Get the time field containing the time values associated
# with the data in the time-enabled layer
startTimeField = str(lyrTime.startTimeField)

# Specify the output mosaic dataset to which the interpolated rasters will be added 
outputMosaicDataset = r"C:\Users\Paul\Documents\ArcGIS\Python Time\Time Data\Time_Enabled.shp"

i = 0
while i <= iterations:
    # Formulate the time query and increment the time by the timeStepInterval
    currentTime = str(startTime + (i*timeStepInterval))
    timeQuery = "\"" + startTimeField + "\"" + " = date '" + currentTime + "'"

    # Create an in-memory feature layer containing points that are valid at each timestep
    tempFeatureLyr = "tempTimeLayer" + str(i)
    arcpy.MakeFeatureLayer_management(lyr, tempFeatureLyr, timeQuery)

    # Create an interpolated raster surface using the points valid at each timestep
    outRaster = r"C:\Users\Paul\Documents\ArcGIS\Python Time\Time Data\Time Data.gdb\raster" + str(i)
    print outRaster
    arcpy.sa.Idw(tempFeatureLyr, "Temp", outRaster)

    # Add the newly created raster surface to a Mosaic Dataset
    arcpy.AddRastersToMosaicDataset_management(outputMosaicDataset, "Raster Dataset", outRaster)

    i = i + 1
    
# Calculate the statistics on the output Mosaic Dataset for
# classifying your data after new rasters are added   
arcpy.CalculateStatistics_management(outputMosaicDataset,"1","1","#")

MathewCoyle · ‎02-04-2013

Are you sure your layer is time enabled?

PaulAlexander1 · ‎02-04-2013

I really think so - at least time slider seems to work OK in ArcMap. Originally was trying to use hourly data, but there is an issue with Ireland to US dates, so below is hourly data but I converted to dates (so from the 1st to the 2nd of a given month is actually from hour 1 to hour 2 - conceptually, I don't need to preserve the actual times, so long as I end up with 24 IDWs)

From the attributes:

OBJECTID Object ID
Station Text
startTime Date
endTime Date
Temp Double
latitute Double
longitude Double

e.g.

FID Shape * OBJECTID Station startTime endTime Temp latitute longitude
0 Point 1 GPEP #1 01/01/2011 02/01/2011 8.9 53.306238 -6.222935
1 Point 2 GPEP #1 02/01/2011 03/01/2011 8.7 53.306238 -6.222935
2 Point 3 GPEP #1 03/01/2011 04/01/2011 8.5 53.306238 -6.222935
3 Point 4 GPEP #1 04/01/2011 05/01/2011 8.1 53.306238 -6.222935
4 Point 5 GPEP #1 05/01/2011 06/01/2011 9.2 53.306238 -6.222935
5 Point 6 GPEP #1 06/01/2011 07/01/2011 10.8 53.306238 -6.222935
6 Point 7 GPEP #1 07/01/2011 08/01/2011 12.2 53.306238 -6.222935
7 Point 8 GPEP #1 08/01/2011 09/01/2011 12.8 53.306238 -6.222935
8 Point 9 GPEP #1 09/01/2011 10/01/2011 14.1 53.306238 -6.222935
9 Point 10 GPEP #1 10/01/2011 11/01/2011 16.7 53.306238 -6.222935
10 Point 11 GPEP #1 11/01/2011 12/01/2011 17.3 53.306238 -6.222935
11 Point 12 GPEP #1 12/01/2011 13/01/2011 17.7 53.306238 -6.222935
12 Point 13 GPEP #1 13/01/2011 14/01/2011 17.7 53.306238 -6.222935
13 Point 14 GPEP #1 14/01/2011 15/01/2011 17.9 53.306238 -6.222935
14 Point 15 GPEP #1 15/01/2011 16/01/2011 19.1 53.306238 -6.222935
15 Point 16 GPEP #1 16/01/2011 17/01/2011 19.1 53.306238 -6.222935
16 Point 17 GPEP #1 17/01/2011 18/01/2011 18.8 53.306238 -6.222935
17 Point 18 GPEP #1 18/01/2011 19/01/2011 18.1 53.306238 -6.222935
18 Point 19 GPEP #1 19/01/2011 20/01/2011 18.7 53.306238 -6.222935
19 Point 20 GPEP #1 20/01/2011 21/01/2011 17.1 53.306238 -6.222935
20 Point 21 GPEP #1 21/01/2011 22/01/2011 17.1 53.306238 -6.222935
21 Point 22 GPEP #1 22/01/2011 23/01/2011 18 53.306238 -6.222935
22 Point 23 GPEP #1 23/01/2011 24/01/2011 18 53.306238 -6.222935
23 Point 24 GPEP #1 24/01/2011 01/01/2011 17.3 53.306238 -6.222935
24 Point 25 GPEP #2 01/01/2011 02/01/2011 10.8 53.292258 -6.138301
25 Point 26 GPEP #2 02/01/2011 03/01/2011 10.6 53.292258 -6.138301
26 Point 27 GPEP #2 03/01/2011 04/01/2011 10.2 53.292258 -6.138301
27 Point 28 GPEP #2 04/01/2011 05/01/2011 10.1 53.292258 -6.138301
28 Point 29 GPEP #2 05/01/2011 06/01/2011 10.5 53.292258 -6.138301
29 Point 30 GPEP #2 06/01/2011 07/01/2011 11.4 53.292258 -6.138301
30 Point 31 GPEP #2 07/01/2011 08/01/2011 11.8 53.292258 -6.138301
31 Point 32 GPEP #2 08/01/2011 09/01/2011 12.4 53.292258 -6.138301
32 Point 33 GPEP #2 09/01/2011 10/01/2011 12.9 53.292258 -6.138301
33 Point 34 GPEP #2 10/01/2011 11/01/2011 14.5 53.292258 -6.138301
34 Point 35 GPEP #2 11/01/2011 12/01/2011 16.4 53.292258 -6.138301
35 Point 36 GPEP #2 12/01/2011 13/01/2011 17.2 53.292258 -6.138301
36 Point 37 GPEP #2 13/01/2011 14/01/2011 16.7 53.292258 -6.138301
37 Point 38 GPEP #2 14/01/2011 15/01/2011 17.3 53.292258 -6.138301
38 Point 39 GPEP #2 15/01/2011 16/01/2011 16.8 53.292258 -6.138301
39 Point 40 GPEP #2 16/01/2011 17/01/2011 17.1 53.292258 -6.138301
40 Point 41 GPEP #2 17/01/2011 18/01/2011 17.1 53.292258 -6.138301
41 Point 42 GPEP #2 18/01/2011 19/01/2011 16.8 53.292258 -6.138301

PaulAlexander1 · ‎02-04-2013

I know you can insert arguments into the script so that users have to define fields, layers etc. when they run the script - could that be a way to work around the problem or is it very complex?

Zeke · ‎02-04-2013

Are you sure you have time enabled, not just time fields in your data (see image below). If time slider is working, I'd think you do, but worth a check.
Still pretty sure the arcpy.gp.IDW_sa call is incorrect, even if it's from help. There are typos there too. Check out the IDW samples in help, they follow Matthew's format Well, they import the whole sa module, but same idea. I'd link, but getting a proxy error at the site right now.

User entered parameters are generally entered when you add a script to a toolbox. Not hard, but I don't think that's the issue here. If it won't run with stuff hard coded, it won't run with parameters either.

PaulAlexander1 · ‎02-04-2013

Thanks for your reply,

It is Definitely time enabled as per your screen-grab there.

Could anyone suggest an alternative (less tidy) approach? Instead of using time stamps a more basic batch running of IDW based on a defined variable (Z i.e. the basis for the IDW) and a matching code (hour one could be coded 1, hour two coded 2) and so forth?

- Paul.