烦死了老报错！！！

最新推荐文章于 2024-08-28 16:14:15 发布

GUAI_GUAI_daydayup

最新推荐文章于 2024-08-28 16:14:15 发布

阅读量81

点赞数

文章标签： python numpy 机器学习

本文链接：https://blog.csdn.net/qq_58291072/article/details/127521216

版权



import os
import re  # regular expressions
import warnings
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
import rasterio as rio
from rasterio.plot import plotting_extent
from osgeo import gdal, osr
print('///')
# os.environ['PROJ_LIB'] =r'C:\Users\jingjing\anaconda3\Library\share\proj'
os.environ['PROJ_LIB'] = 'D:\ProgramData\Anaconda3\Lib\site-packages\pyproj\proj_dir\share\proj'
def GetFullFileName(filePath):
    fileList = os.listdir(filePath)
    fullList = []
    for file in fileList:
        fullList.append(os.path.join(filePath,file))
    return fileList,fullList
print('llllllllllllllllllll')
def OpenHDF(fileName, Var_name):
    with rio.open(fileName) as dataset:
        for name in dataset.subdatasets:
            #Use regular expression to identify if subdataset has LST_Day_1km in the name (the bands)
            if re.search(Var_name,name):
                with rio.open(name) as subdataset:
                    mod = subdataset.read(1)
                    mod_meta = subdataset.meta
    return mod,mod_meta
print('======================== ')
def CreateTiff(fullFileName,dynamicArray, width, length, proj,geotransform):
    driver = gdal.GetDriverByName('GTiff')
    out_tif = driver.Create(fullFileName, width, length, 1, gdal.GDT_Float32)  # 创建框架
    out_tif.SetGeoTransform((geotransform[2],geotransform[0],geotransform[1],geotransform[5],geotransform[3],geotransform[4]))
    out_tif.SetProjection(str(proj))  # 给新建图层赋予投影信息
    
    # 数据写出
    out_tif.GetRasterBand(1).WriteArray(dynamicArray)  # 将数据写入内存，此时没有写入硬盘
    out_tif.FlushCache()  # 将数据写入硬盘
    out_tif = None  # 注意必须关闭tif文件

if __name__ == "__main__":
    
    filePath = r'D:\MODIS666\batchmodis_try'
    output_path = r'D:\MODIS666\666'
    
    file_list,full_file_name_list = GetFullFileName(filePath)
    
    for file in file_list:
        Var_name = "1 km 16 days NDVI"
        lst_night,lst_night_meta=OpenHDF(os.path.join(filePath,file), Var_name)
        lst_night = lst_night * 0.02 # NDVI的转换系数为0.0001
        width = (lst_night.shape)[1]
        length = (lst_night.shape)[0]
        
        # print(lc_meta)
        
        fullFileName = os.path.join(output_path,file[0:-4]+'1 km 16 days NDVI.tif')
        print('qqqqqqqqqqqqqqqqqqqqqq')
        CreateTiff(fullFileName,lst_night, width, length, lst_night_meta['crs'], lst_night_meta['transform'])
        
        # print(file)
        # fullFileName1 = os.path.join(output_path,file[0:-4]+'_LST_night_WGS.tif')
        # os.system('gdalwarp -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs" -t_srs "+proj=longlat +datum=WGS84 +no_defs" -srcnodata 0 -dstnodata 0 {} {}'.format(fullFileName, fullFileName1))