Select to view content in your preferred language

Extract NetCDF raster values to 'timestepped' shapefiles

1657
3
10-15-2013 03:33 AM
ChristopherVos1
Emerging Contributor
I have a NetCDF raster for Sea Surface Temperature (SST), which has a time dimension that is monthly from 1950 - 2012.
This is a single NetCDF file with many 'time slices' - 744 i think (12 per year for 62 years)

I have hurricane track data which is in point shape files, which I have also been able to time enable to this same monthly increments.

ArcMap understands this - and I can run animations which has the NetCDF raster and shapefiles changing correctly with time.

What I am trying to do is to now is produce a model (workflow) in model builder which extracts the SST value for each data point at the correct time.

Please provide assistance.

Chris
0 Kudos
3 Replies
ChrisFox3
Frequent Contributor
To accomplish this you are going to need to use python, which you could then create a script tool from and use in model builder. Below is a code sample. Basically we are opening the MXD, gettting some references to the data frame and layers and then stepping through time. At each step we perform our analysis and then move to the next time step.

#Get a reference to the map document and first data frame
mxd = arcpy.mapping.MapDocument('Current') #If code is run in the app with the MXD open, otherwise pass the path to the MXD
df = arcpy.mapping.ListDataFrames(mxd)[0] #Assuming here there is only one data frame in the MXD

rasterLayerName = "SST" #replace with name of netCDF raster layer
rasterLayer = arcpy.mapping.ListLayers(mxd, rasterLayerName, df)[0]

hurricaneLayer = "Hurricanes" #replace with name of hurricane layer
hurricaneLayer = arcpy.mapping.ListLayers(mxd, hurricaneLayerName, df)[0]

#Loop through each time step within the full time extent
df.time.currentTime = df.time.startTime
while df.time.currentTime <= df.time.endTime:

    #Add your logic here, I assume you are using extract values to point
    #So we don't overwrite the output each time, we need a unique name, this will append the Year and Month value currently being processed.
    uniqueName = "Hurricanes_" + dt.time.currentTime.strftime("%Y_%m") 
    outPointFeatures = r"C:\OutputDate\" + uniqueName
    ExtractValuesToPoints(hurricaneLayer, rasterLayer, outPointFeatures,
                      "NONE", "VALUE_ONLY")

    #Set the current time of the data frame to the next time step
    df.time.currentTime += df.time.timeStepInterval
0 Kudos
ChristopherVos1
Emerging Contributor
Chris, Thanks for your reply.

Everything was working in the Python environment until I typed in the line:

hurricaneLayer = arcpy.mapping.ListLayers(mxd, hurricaneLayerName, df)[0]

I received a Runtime Error saying: <type 'exceptions.NameError'> name 'hurricaneLayerName' is not defined

I tried changing the code to:

hurricaneLayer = arcpy.mapping.ListLayers(mxd, hurr_data, df)[0] (hurr_data is the name of the shapefile feature class)

What do i need to do differently?

Kind Regards,

Chris
0 Kudos
ChrisFox3
Frequent Contributor
Sorry there was a typo in the code. The variable above should be hurricaneLayerName.

hurricaneLayerName = "Hurricanes" #replace with name of hurricane layer
hurricaneLayer = arcpy.mapping.ListLayers(mxd, hurricaneLayerName, df)[0]


You were on the right track though. You can just pass in the name, but it needs to be the name of the layer in the map, not the shapefile name (They could be different). You would also need to surrond it in quotes, because it is a string.

hurricaneLayer = arcpy.mapping.ListLayers(mxd, "hurr_data", df)[0]
0 Kudos