自己看，以防遗忘。

记录硕士期间的论文程序代码——1
python

# from mpl_toolkits.basemap import Basemap
from osgeo import gdal, osr, ogr
# import matplotlib.pyplot as plt
import numpy as np
import os
from pyhdf.SD import *
import xlwt
import xlsxwriter
import pandas as pd
# import itertools
import time
from libtiff import TIFF
# a=os.listdir(os.getcwd())
working_path = "D:\\mod\\2020\\h24v05\\"
os.chdir(working_path)
input_path = "D:\\mod\\2020\\h24v05\\"
file_list = os.listdir(input_path)
class SNOWPixcel():
    """docstring for SNOWPixcel"""
    def __init__(self, value):
        self.value = value
        pass

    def GetSnowCoverDuration(self):
        return np.sum(self.value == 10)

    def GetSnowFirstAndLastDay(self):
        first_day, last_day = 0, 0
        indexlists= np.where(self.value == 10)
        index_list=indexlists[0]+1
        gap =0
        max_scd, min_sed = 0, 0
        for i in range(1, len(index_list)):
            if index_list[i] - index_list[i - 1] > gap:
                gap=index_list[i] - index_list[i - 1]
                max_scd, min_sed = index_list[i], index_list[i-1]
#                if min_sed>114:
#                   min_sed+=1
#                if max_scd>114:
 #                   max_scd+=1
#else:
#                    if max_scd>235:
#                        max_scd+=2
 #                   else:
#                        pass
            else:
                pass
        return max_scd, min_sed
col =2400
row = 2400
pixel = col * row
# 
snow_cover = np.zeros((len(file_list), pixel), dtype = np.int)
# 
print('start')
print(time.asctime(time.localtime(time.time())))
# print("start-read" + datetime.now())
for i in range(len(file_list)):
    print("dataset:" + str(i))
    dataset = SD(file_list[i])
    data = dataset.select('NDSI_Snow_Cover').get()
    if data.shape[0] == row and data.shape[1] == col:
        snow_cover[i, :] = data.flatten()
    else:
        print(file_list[i] + 'Wrong Shape')
    del dataset
print('end')
print(time.asctime(time.localtime(time.time())))
# 
snow_cover_mask = snow_cover.copy()
snow_cover[snow_cover_mask == 0] = 500
snow_cover[snow_cover <=100] = 10
# snow cover duration
result = np.zeros((3, pixel), dtype = np.int)
# 
for j in range(snow_cover.shape[1]): # 按列遍历
    print('process' + str(j))
    sp = SNOWPixcel(snow_cover[:, j])
    '''
    snow_time_series = snow_cover[:, j]
    SCD[j] = np.sum(snow_time_series == 10)
    '''
    result[0, j] = sp.GetSnowCoverDuration()
    days = sp.GetSnowFirstAndLastDay()
    result[1, j] = int(days[0])
    result[2, j] = int(days[1])
    pass
B=np.array(result[1,]).reshape(col,row)
SOD_reshape=pd.DataFrame(B)
writer=pd.ExcelWriter("D:\\mod\\2020\\2020h24v05SOD.xlsx")		# 写入Excel文件
SOD_reshape.to_excel(writer, 'page_1', float_format='%.5f')		# page_1是写入excel的sheet名
writer.close()
C=np.array(result[0, ]).reshape(col,row)
SCD_reshape=pd.DataFrame(C)
writer=pd.ExcelWriter("D:\\mod\\2020\\2020h24v05SCD.xlsx")		# 写入Excel文件
SCD_reshape.to_excel(writer, 'page_1', float_format='%.5f')		# page_1是写入excel的sheet名
writer.close()
D=np.array(result[2, ]).reshape(col,row)
SED_reshape=pd.DataFrame(D)
writer=pd.ExcelWriter("D:\\mod\\2020\\2020h24v05SED.xlsx")		# 写入Excel文件
SED_reshape.to_excel(writer, 'page_1', float_format='%.5f')		# page_1是写入excel的sheet名
writer.close()
tif = TIFF.open("D:\\mod\\2020\\2020h24v05SOD.tif", 'w') 
tif.write_image(SOD_reshape.astype(np.uint16), compression=None)
tif.close() #--> 16bit(0~65535,8位则0~255)
tif = TIFF.open("D:\\mod\\2020\\2020h24v05SCD.tif", 'w') 
tif.write_image(SCD_reshape.astype(np.uint16), compression=None)
tif.close() #--> 16bit(0~65535,8位则0~255)
tif = TIFF.open("D:\\mod\\2020\\2020h24v05SED.tif", 'w') 
tif.write_image(SED_reshape.astype(np.uint16), compression=None)
tif.close() #--> 16bit(0~65535,8位则0~255)