I have written code using Python for Fully Constrained Least Squares (FCLS) Linear Spectral Mixture Analysis, which could be applied for unmixing multispectral image successfully.
However, the operation efficiency is very low, taking about 2 hours for each MODIS image (rows: 1620, columns: 3024, bands：7, blocksize: 512*512). Is there any method for optimization this code to improve efficiency?
# -*- coding: utf-8 -*-
# author = qiangsun
import pysptools.abundance_maps as amp
import arcpy
from arcpy import env
import numpy as np
import os
import datetime
# path of modis reflectance images
ref_path = ["I:\MODIS_Processing\CS\CE"]
for image in ref_path:
# workspace
env.workspace = image
LIST = arcpy.ListRasters("*", "dat")
# output result
fileout = os.path.join("I:\MODIS_Processing\sma_block",r"output/block1.tif")
# Set environment for output
arcpy.env.overwriteOutput = True
# size
blocksize = 512
# EMs
GV =[0.023978245,0.076136685,0.05170522,0.479961538,0.453641505,0.251023205,0.101172368]
DA = [0.048255164,0.069696024,0.0850382,0.10709865,0.1255218,0.12941472,0.130492952]
SL = [0.097499328,0.174871013,0.240758917,0.29218451,0.360858768,0.403369335,0.424446828]
SA = [0.256721226,0.373879832,0.44604235,0.500752414,0.436125654,0.388239232,0.205382038]
IC = [0.923212177,0.940795633,0.964715637,0.981163445,0.577684378,0.15711141,0.069484282]
EM = np.array([GV,SL,SA,IC,DA])
EM_mat = np.mat(EM)
for myRaster1 in LIST:
arcpy.env.outputCoordinateSystem = myRaster1
starttime = datetime.datetime.now()
print "running sma of %s" % myRaster1
myRaster = arcpy.Raster(myRaster1.encode('ascii'))
filelist1 = []
filelist2 = []
filelist3 = []
filelist4 = []
filelist6 = []
filelist7 = []
blockno1 = 0
blockno2 = 0
blockno3 = 0
blockno4 = 0
blockno6 = 0
blockno7 = 0
fileout1 = "I:\MODIS_Processing\sma_block\SMA_%s_GV.tif" % myRaster1[:-4]
fileout2 = "I:\MODIS_Processing\sma_block\SMA_%s_DA.tif" % myRaster1[:-4]
fileout3 = "I:\MODIS_Processing\sma_block\SMA_%s_SL.tif" % myRaster1[:-4]
fileout4 = "I:\MODIS_Processing\sma_block\SMA_%s_SA.tif" % myRaster1[:-4]
fileout6 = "I:\MODIS_Processing\sma_block\SMA_%s_IC.tif" % myRaster1[:-4]
fileout7 = "I:\MODIS_Processing\sma_block\SMA_%s_RMS.tif" % myRaster1[:-4]
for x in range(0, myRaster.width, blocksize):
for y in range(0, myRaster.height, blocksize):
mx = myRaster.extent.XMin + x * myRaster.meanCellWidth
my = myRaster.extent.YMin + y * myRaster.meanCellHeight
# Upper right coordinate of block (in cells)
nx = min([x + blocksize, myRaster.width])
ny = min([y + blocksize, myRaster.height])
B = arcpy.RasterToNumPyArray(myRaster, arcpy.Point(mx, my), ncols=nx - x, nrows=ny - y)
GVA = []
DAA = []
SLA = []
SAA = []
ICA = []
RMSA = []
for i in range(0, B.shape[1], 1):
for j in range(0, B.shape[2], 1):
b = B[:, i, j]
b1 = []
for m in range(0, len(b), 1):
mn = float(b[m])
b1.append(mn)
var_y = np.mat(b1)
# processing of SMA
SMA = amp.amaps.FCLS(var_y, EM_mat)
FR_GV = SMA[0][0]
FR_SL = SMA[0][1]
FR_SA = SMA[0][2]
FR_IC = SMA[0][3]
FR_DA = SMA[0][4]
f0 = (b1[0] - (FR_GV * GV[0] + FR_DA * DA[0] + FR_SL * SL[0] + FR_SA * SA[0] + FR_IC * IS[0])) ** 2
f1 = (b1[1] - (FR_GV * GV[1] + FR_DA * DA[1] + FR_SL * SL[1] + FR_SA * SA[1] + FR_IC * IS[1])) ** 2
f2 = (b1[2] - (FR_GV * GV[2] + FR_DA * DA[2] + FR_SL * SL[2] + FR_SA * SA[2] + FR_IC * IS[2])) ** 2
f3 = (b1[3] - (FR_GV * GV[3] + FR_DA * DA[3] + FR_SL * SL[3] + FR_SA * SA[3] + FR_IC * IS[3])) ** 2
f4 = (b1[4] - (FR_GV * GV[4] + FR_DA * DA[4] + FR_SL * SL[4] + FR_SA * SA[4] + FR_IC * IS[4])) ** 2
f5 = (b1[5] - (FR_GV * GV[5] + FR_DA * DA[5] + FR_SL * SL[5] + FR_SA * SA[5] + FR_IC * IS[5])) ** 2
f6 = (b1[6] - (FR_GV * GV[6] + FR_DA * DA[6] + FR_SL * SL[6] + FR_SA * SA[6] + FR_IC * IS[6])) ** 2
RMS = np.sqrt((f1 + f2 + f3 + f4 + f5 + f6 + f0) / 8)
GVA.append(FR_GV)
DAA.append(FR_DA)
SLA.append(FR_SL)
SAA.append(FR_SA)
ICA.append(FR_IC)
RMSA.append(RMS)
gv_a = np.array(GVA).reshape(B.shape[1], B.shape[2])
da_a = np.array(DAA).reshape(B.shape[1], B.shape[2])
sl_a = np.array(SLA).reshape(B.shape[1], B.shape[2])
sa_a = np.array(SAA).reshape(B.shape[1], B.shape[2])
ic_a = np.array(ICA).reshape(B.shape[1], B.shape[2])
rms_a = np.array(RMSA).reshape(B.shape[1], B.shape[2])
GV_A = arcpy.NumPyArrayToRaster(gv_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
DA_A = arcpy.NumPyArrayToRaster(da_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
SL_A = arcpy.NumPyArrayToRaster(sl_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
SA_A = arcpy.NumPyArrayToRaster(sa_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
IC_A = arcpy.NumPyArrayToRaster(ic_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
RMS_A = arcpy.NumPyArrayToRaster(rms_a, arcpy.Point(mx, my), myRaster.meanCellWidth,
myRaster.meanCellHeight)
filetemp1 = ('1_%i.' % blockno1).join(fileout1.rsplit('.', 1))
filetemp2 = ('2_%i.' % blockno2).join(fileout2.rsplit('.', 1))
filetemp3 = ('3_%i.' % blockno3).join(fileout3.rsplit('.', 1))
filetemp4 = ('4_%i.' % blockno4).join(fileout4.rsplit('.', 1))
filetemp6 = ('6_%i.' % blockno6).join(fileout6.rsplit('.', 1))
filetemp7 = ('7_%i.' % blockno7).join(fileout7.rsplit('.', 1))
GV_A.save(filetemp1)
DA_A.save(filetemp2)
SL_A.save(filetemp3)
SA_A.save(filetemp4)
IC_A.save(filetemp6)
RMS_A.save(filetemp7)
filelist1.append(filetemp1)
blockno1 += 1
filelist2.append(filetemp2)
blockno2 += 1
filelist3.append(filetemp3)
blockno3 += 1
filelist4.append(filetemp4)
blockno4 += 1
filelist6.append(filetemp6)
blockno6 += 1
filelist7.append(filetemp7)
blockno7 += 1
filename1 = "%s_GV.tif" % myRaster1[:-4]
filename2 = "%s_DA.tif" % myRaster1[:-4]
filename3 = "%s_SL.tif" % myRaster1[:-4]
filename4 = "%s_SA.tif" % myRaster1[:-4]
filename6 = "%s_IS.tif" % myRaster1[:-4]
filename7 = "%s_RMSE.tif" % myRaster1[:-4]
arcpy.MosaicToNewRaster_management(";".join(filelist1[0:]), "I:/MODIS_Processing/SMA/", filename1, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
arcpy.MosaicToNewRaster_management(";".join(filelist2[0:]), "I:/MODIS_Processing/SMA/", filename2, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
arcpy.MosaicToNewRaster_management(";".join(filelist3[0:]), "I:/MODIS_Processing/SMA/", filename3, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
arcpy.MosaicToNewRaster_management(";".join(filelist4[0:]), "I:/MODIS_Processing/SMA/", filename4, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
arcpy.MosaicToNewRaster_management(";".join(filelist6[0:]), "I:/MODIS_Processing/SMA/", filename6, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
arcpy.MosaicToNewRaster_management(";".join(filelist7[0:]), "I:/MODIS_Processing/SMA/", filename7, " ", "64_BIT",
"0.00417262",
"1", "FIRST", "FIRST")
for fileitem1 in filelist1:
if arcpy.Exists(fileitem1):
arcpy.Delete_management(fileitem1)
for fileitem2 in filelist2:
if arcpy.Exists(fileitem2):
arcpy.Delete_management(fileitem2)
for fileitem3 in filelist3:
if arcpy.Exists(fileitem3):
arcpy.Delete_management(fileitem3)
for fileitem4 in filelist4:
if arcpy.Exists(fileitem4):
arcpy.Delete_management(fileitem4)
for fileitem6 in filelist6:
if arcpy.Exists(fileitem6):
arcpy.Delete_management(fileitem6)
for fileitem7 in filelist7:
if arcpy.Exists(fileitem7):
arcpy.Delete_management(fileitem7)
del myRaster
endtime = datetime.datetime.now()
print "%s:%s"%(myRaster1,(endtime - starttime).seconds)
print "run over"
amp.amaps.FCLS(var_y, EM_mat)but it's an educated guess at best. Running the code under a profiler, such as cProfile, is probably the way to go here. Nothing obvious to me concerning performance improvement, anyway. \$\endgroup\$