- After having converted the original monthly MOD13C2 product (from year 2000-2020) to raster format, now I have to calculate the mean value of the 251 rasters.
First I have tried this tutorial (https://blog.csdn.net/lucky51222/article/details/109412616) which simply used the algebra function that behaved badly in cooperating NA values. So I tried the second one (https://blog.csdn.net/qq_35056050/article/details/112017543) which converted each raster to array to skip NA pixels. I adapted it into my code:
import arcpy, sys, os, glob from arcpy.sa import * import numpy arcpy.CheckOutExtension('Spatial') # input path inws = "G:/data0610/MODIS_VI/EVI/EVI_pro/" # output path outws = "G:/data0610/MODIS_VI/EVI/EVI_pro/" rasters = glob.glob(os.path.join(inws, "*.tif")) r = Raster(rasters[0]) array = arcpy.RasterToNumPyArray(r) # convert to numpy rowNum, colNum = array.shape sum = numpy.zeros(shape=array.shape) # save the accumulating value count = numpy.zeros(shape=array.shape) # save the counting number Average = numpy.zeros(shape=array.shape) # save the mean value for ras in rasters: rmm = Raster(ras) array = arcpy.RasterToNumPyArray(rmm) # one by one pixel for i in range(0, rowNum): for j in range(0, colNum): if array[i][j] >= 0 : # verdict invalid value sum[i][j] += array[i][j] # accumulate count[i][j] += 1 # counter continue Average = sum / count # cal the mean value # save the raster lowerLeft = arcpy.Point(r.extent.XMin, r.extent.YMin) cellWidth = r.meanCellWidth cellHeight = r.meanCellHeight nameT = "evi.tif" outname = os.path.join(outws, nameT) arcpy.env.overwriteOutput = True #convert to WGS84 inf = "G:/data0610/MODIS_VI/sm_mean.tif" arcpy.env.outputCoordinateSystem = Raster(inf) # convert the crs to wgs84 print("successfully converted the CRS!") AvgRas = arcpy.NumPyArrayToRaster(Average, lowerLeft, cellWidth, cellHeight, r.noDataValue) # turn into raster AvgRas.save(outname) print("successfully output the evi_mean.tif!") # resample outname_res = outws + "evi_mean_res.tif" # get the standard cellsize cellsize025 = "{0} {1}".format(arcpy.Describe(inf).meanCellWidth, arcpy.Describe(inf).meanCellHeight) arcpy.Resample_management(AvgRas, outname_res, cellsize025, "NEAREST") print("successfully output the evi_mean.tif with the 0.25 degree resolution!")
Unfortunately, the arcpy (python 2.7, 32 bit) had the memory error because there were too many large rasters (I'm dealing with the global extent). I found the reason and solution from this question (https://gis.stackexchange.com/questions/291948/numpy-memory-error-on-large-rasters). (See pic.1)
So I installed the 64-bit Background Geoprocessing of ArcGIS and run the above code again, then it came up with another problem. (See pic.2)
It turned out the soultion might be useless, because ArcGIS is very sensitive to administrator license, you can't run the 64-bit python while the ArcGIS being 32-bit.
Now return to my initial question: How to calculate the mean value of multiple rasters in python/arcpy? Did I make things complicated? Is there a simpler way of generating the mean value raster?
It's really driving me crazy.
- Badges
- Science topic
- Similar topics
- Biological Science
- Bioecology
- Systems Ecology
- Landscape Ecology
- Connectivity