2
\$\begingroup\$

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"
\$\endgroup\$
  • 1
    \$\begingroup\$ It's hard to tell where most time is spent especially not knowing well the framework. I believe most time is spent in 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\$ – Arthur Havlicek Nov 17 '18 at 3:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.