Creating a list of raster lists

[JamesKelly2]

Hello,

I am a beginner at python programming, and this is my first post to GeoNet. So here it goes!

I have a huge list of rasters downloaded from the SNODAS website for years 2003 - present (intra-annual range of 11/1-3/31). I have a .tif file for everyday during this period which contains the snow depth data for that day. The metadata for each raster is stored in the tif file's name(e.g., us_ssmv11036tS__T0001TTNATS2013111305HP001.dat.reproj.tif ). However, I am only interested in the date (Nov. 13, 2013 in this case).

My objective is to split this huge list of rasters into winter-years (i.e., rasters from Jan-Mar of a given year belong to the winter-year of the previous calendar year) so that I can pass each winter-year list of rasters through the greaterThanFrequency tool with a simple loop.

Below is a script that works if I download one year of data at a time. I'd appreciate any feedback on how get a list of rasters for each winter-year.

import arcpy
from arcpy import env
from arcpy.sa import *


env.workspace = "E:\DECGIS\snow\NSIDC\projects" \
"\Searchlight\production\data\NSIDC_Data_1021063617170" \
"\Snow_Data_Assimilation_System_(SNODAS)_Data_Products_at_NSIDC"
env.overwriteOutput = True
arcpy.env.snapRaster = "us_ssmv11036tS__T0001TTNATS2013111305HP001.dat.reproj.tif"
arcpy.env.outputCoordinateSystem = "us_ssmv11036tS__T0001TTNATS2013111305HP001.dat.reproj.tif"
arcpy.env.nodata = "-9999"
inRasters = arcpy.ListRasters("*", "TIF")
arcpy.CheckOutExtension("Spatial")

#create raster of constant values for greaterThanFrequency tool
constantValue = "381"
outpath3 = str(env.workspace+"\\threshold4.tif")
cellSize = '8.33333333813986E-03'
outConstRaster = CreateConstantRaster(constantValue, 'INTEGER', cellSize,"us_ssmv11036tS__T0001TTNATS2013111305HP001.dat.reproj.tif")
outConstRaster.save(outpath3)

#run greater than frequency analysis
inValueRaster = "threshold4.tif"
outGTF = GreaterThanFrequency(inValueRaster, inRasters)

# Save the output 
outGTF.save(r"E:\DECGIS\snow\output\2013daysGt15.tif")

[GerryGabrisch]

James,

Assuming you have all the tiffs in one directory, and assuming the tiff file names are of the same length, you could sift through all the files in the folder, find which files were winter months, create a list of those file names, and then execute you code after the list of file names is complete.  This gets the file names, splits and slices out the month, and writes the file name to a new list.

import  os
indir = ""#YOUR DIRECTORY OF TIFFS

#Create an empty list to store all the file names for winter months....
listofdatatoprocess = []

#Get the files in the directory....
for root, directories, files in os.walk(indir):
    for filename in files:
      filenamesplit = filename.split(".")
        #extract the month...    
      month = filenamesplit[0][31:33]
      if month == '11' or month == '12' or month == '01' or month == '02' or month == '03':
            listofdatatoprocess.append(filename)

for item in listofdatatoprocess:
        #YOUR CODE HERE

[NeilAyres]

And you could replace that if month == ...

With

if month in ['11','12','01','02','03'] :

[JamesKelly2]

Hi Gerry and Neil,

Thanks for looking into my issue.

Either I was unclear in my explanation or I am not understanding your explanation.

I am not sure how this will help me make a separate list for each year's winter in the dataset. There are 11 years worth of data in the same folder. The only months I downloaded for each year are 11,12,1,2 and 3 so I don't need to separate out winter months. What I need to separate are the years(e.g., winter of '03-'04). This is a challenge because months 1, 2, and 3 need to be assigned to the previous calendar year.

Any thoughts?

Thanks again!

[curtvprice · Accepted solution]

Here's a simple approach, generate a list of year/month/days, and then convert the list into raster paths. We could get fancier with the datetime module but I think this should work fine. I set this up for 2003 / 2004 but you could nest this inside a loop to process each year at a time.

# us_ssmv11036tS__T0001TTNATS2013111305HP001.dat.reproj.tif
daylist = []
year = 2003
for mon in [11, 12]:
    for day in range(1,32):
        daylist.append("{:04d}{:02d}{:02d}".format(year,mon,day))
year += 1
for mon in [1,2,3]:
    for day in range(1,32):
        daylist.append("{:04d}{:02d}{:02d}".format(year,mon,day))
# convert the day list into tif pathnames
tpath = "us_ssmv11036tS__T0001TTNATS{}05HP001.dat.reproj.tif"   
rpaths = [tpath.format(d) for d in daylist]
# if the raster exists, create a raster object and add to the list
# this raster list is the one you can pass to your tool
rasters = [Raster(rpath) for rpath in rpaths if arcpy.Exists(rpath)]

[JamesKelly2]

Thanks Curtis! This is exactly what I need!