Calculate daily average of raster files based on their names in python

3712
27
Jump to solution
06-13-2017 10:15 AM
AmenehTavakol
New Contributor III

I have hundreds of SMAP soil moisture data and need to calculate the daily average of them. For every day there are 8 files (every 3 hours) and I need a code to calculate the average based on their names. 

I have a code but i am interested to know how can I define a code that average data based on their names? 

# Import arcpy module
import arcpy


# Local variables:
SMAP_L4_SM_gph_20150331T043000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c15_tif = "SMAP_L4_SM_gph_20150331T043000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c15.tif"
SMAP_L4_SM_gph_20150331T073000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0d_tif = "SMAP_L4_SM_gph_20150331T073000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0d.tif"
SMAP_L4_SM_gph_20150331T103000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c14_tif = "SMAP_L4_SM_gph_20150331T103000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c14.tif"
SMAP_L4_SM_gph_20150331T133000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf5_tif = "SMAP_L4_SM_gph_20150331T133000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf5.tif"
SMAP_L4_SM_gph_20150331T163000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf3_tif = "SMAP_L4_SM_gph_20150331T163000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf3.tif"
SMAP_L4_SM_gph_20150331T193000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf4_tif = "SMAP_L4_SM_gph_20150331T193000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf4.tif"
SMAP_L4_SM_gph_20150331T223000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0f_tif = "SMAP_L4_SM_gph_20150331T223000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0f.tif"
SMAP_L4_SM_gph_20150331T013000_Vv2030_001_Geophysical_Data_sm_surface_bdea9eb8_tif = "SMAP_L4_SM_gph_20150331T013000_Vv2030_001_Geophysical_Data_sm_surface_bdea9eb8.tif"
v_name_ = "D:\\NASA\\Data\\SMAP\\2014\\Processed\\%name%.tif"

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("( \"%SMAP_L4_SM_gph_20150331T043000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c15.tif%\" + \"%SMAP_L4_SM_gph_20150331T073000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0d.tif%\" + \"%SMAP_L4_SM_gph_20150331T103000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c14.tif%\" +\"%SMAP_L4_SM_gph_20150331T133000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf5.tif%\" + \"%SMAP_L4_SM_gph_20150331T163000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf3.tif%\" + \"%SMAP_L4_SM_gph_20150331T193000_Vv2030_001_Geophysical_Data_sm_surface_bdea9bf4.tif%\" + \"%SMAP_L4_SM_gph_20150331T223000_Vv2030_001_Geophysical_Data_sm_surface_bdea9c0f.tif%\" + \"%SMAP_L4_SM_gph_20150331T013000_Vv2030_001_Geophysical_Data_sm_surface_bdea9eb8.tif%\") / 8", v_name_)


Your help is really appreciated.
0 Kudos
27 Replies
AmenehTavakol
New Contributor III

I put the final codes to help others with the same question:

import arcpy
from shutil import copy2
from collections import defaultdict
arcpy.CheckOutExtension("Spatial")

def main():
    raster_dir = r"D:\NASA\Data\SMAP\OutPut"
    arcpy.env.workspace = raster_dir
    raster_days = defaultdict(list)
    for r in arcpy.ListRasters():
        raster_days[r[15:23]].append(r)
    for date, rasters in raster_days.items():
        print(date)
        outCellStatistics = arcpy.sa.CellStatistics(rasters, "MEAN")
        copy2(r"D:\NASA\Data\SMAP\OutPut\meanc_ras.tif",
              r"D:\NASA\Data\SMAP\OutPut\Processed\meanc_ras_{}.tif".format(date))

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
BlakeTerhune
MVP Regular Contributor

In the future, try Posting code with Syntax Highlighting on GeoNet. Makes code easier to read and copy.

0 Kudos
curtvprice
MVP Esteemed Contributor

I went ahead and fixed it.

0 Kudos
AmenehTavakol
New Contributor III

Dear all, 

This is my last code: I need to copy my output file in a folder:

import arcpy
from shutil import copy2

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

from collections import defaultdict

def main():
    raster_dir = r"D:\NASA\Data\SMAP\OutPut"
    arcpy.env.workspace = raster_dir
    raster_days = defaultdict(list)
    for r in arcpy.ListRasters():
        raster_days[r[15:23]].append(r)

    for date, rasters in raster_days.items():
        print(date)
        #for r in rasters:
            #print("\t{}".format(r))
        outCellStatistics = arcpy.sa.CellStatistics(rasters, "MEAN")
        #outname = r"D:\NASA\Data\SMAP\OutPut\Processed"
        #outCellStatistics.save(outname)
        copy2("D:\NASA\Data\SMAP\OutPut\meanc_ras.tif","D:\NASA\Data\SMAP\OutPut\Processed\meanc_ras_{}.tif".format(date))
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

But now I receive an error :


IOError: [Errno 2] No such file or directory: 'D:\\NASA\\Data\\SMAP\\OutPut\\meanc_ras.tif'

I don't know how to save my results!

0 Kudos
BlakeTerhune
MVP Regular Contributor

Looking at the code samples, try

outCellStatistics.save(r"D:\NASA\Data\SMAP\OutPut\Processed\outCellStatsImage.tif")‍‍‍‍‍
LarryBiodun
New Contributor III

Hello everyone.  I am trying to do a similar thing to the original post however, I have thousands of daily files that I want to use the cell statistics to calculate the 8 day averages. Simply put instead of the final code above doing a mean of the daily values, I want the code to do a mean after every eight days and save each file with the beginning date of each 8 day average. Your help is greatly appreciated!

0 Kudos
XanderBakker
Esri Esteemed Contributor

Hi Larry Biodun , it is better to create a new question and provide more information on the input data (names for instance and storage type) and if possible attach some data. Depending on the names we can see if a similar approach can be applied to your problem.

When you create a new question do include a link to this thread.

0 Kudos
LarryBiodun
New Contributor III

Thanks Xander. I'll do that now

0 Kudos