getVIref.py 文件源码

def import_all_year_data(tile):
    temp = np.zeros([46, 2400*2400], np.dtype(float))
    if int(tile[5:6])<2:
        temp[:]=np.nan
    for doy in range(1, 369, 8):
        evifile = buildVrtFile(root, doy, tile, 'evi')
        cloudfile = buildVrtFile(root, doy, tile, 'cloudmask')
        aerosolfile = buildVrtFile(root, doy, tile, 'aerosolmask')
        #if no file found for this DOY
        if evifile == 0: continue
        #doyList.append(doy)
        #build vrt for EVI
        vrtEVI = os.path.join(os.path.dirname(evifile), str(1000+doy)[1:]+tile+'EVI_vrt.vrt')
        print "Building the vrt file: ", evifile
        os.system('gdalbuildvrt -separate -input_file_list '+evifile+' '+vrtEVI)
        inEVI = gdal.Open(vrtEVI)
        EVI = inEVI.ReadAsArray()
        #build vrt for cloudmask
        vrtcloud = os.path.join(os.path.dirname(cloudfile), str(1000+doy)[1:]+tile+'cloud_vrt.vrt')
        print "Building the vrt file: ", cloudfile
        os.system('gdalbuildvrt -separate -input_file_list '+cloudfile+' '+vrtcloud)
        incloud = gdal.Open(vrtcloud)
        cloud = incloud.ReadAsArray()
        #build vrt for aerosol
        vrtaerosol = os.path.join(os.path.dirname(aerosolfile), str(1000+doy)[1:]+tile+'aerosol_vrt.vrt')
        print "Building the vrt file: ", aerosolfile
        os.system('gdalbuildvrt -separate -input_file_list '+aerosolfile+' '+vrtaerosol)
        inaerosol = gdal.Open(vrtaerosol)
        aerosol = inaerosol.ReadAsArray()
        global rows, cols, geoProj, geoTran
        rows = 2400
        cols = 2400
        geoTran = inEVI.GetGeoTransform()
        geoProj = inEVI.GetProjection()
        #mask for bad quality
        EVIgood = ma.masked_where((cloud != 1)|(aerosol == 0)|(EVI < 0)|(EVI > 10000), EVI)
        EVIgood = EVIgood.reshape(EVIgood.size/2400/2400, 2400*2400)
        medianEVI = np.nanmedian(EVIgood, axis=0)
        EVI = None
        aerosol = None
        cloud = None
        EVIgood = None
        #assign to the 46 layer of matrix
        temp[(doy-1)/8, :] = medianEVI
        meanEVI = None
    return temp