CIA和LGC数据处理成tif并显示的python程序_cia lgc

作者：从前慢现在也慢 | 2024-03-27 08:22:22

踩

cia lgc

CIA和LGC数据是做数据融合的实验数据，一个是modis，一个是landsat，以下是根据不同数据做的python程序，我也是经常在网上参考大家的代码，包括以下的代码，也是参考了网上的信息，希望能给大家带来一点帮助，祝大家学习快乐！
以下是LGC_int转tif的程序


# -*- coding: utf-8 -*-
"""
Created on Thu Sep 28 14:47:13 2023
@author: Administrator
"""
 
# -*- coding: utf-8 -*-
import numpy as np
from osgeo import gdal
import os
 
 
# 依据BSQ存储规则，按照存储完单波段整幅图像后再存储下一波段的存储方法进行提取并存入数组。
def read_as_bsq(dataarr,bands, rows, col):
    imgarr = np.zeros((bands,rows,col))
    for b in range(bands):              #取出每个波段
        start = b * rows * col
        end = start + rows * col
        arrtem = dataarr[start:end]
        for r in range(rows):           #一维数组按行取出，存入对应三维数组。
            start2 = r * col
            end2 = start2 + col
            imgarr[b,r,:] = arrtem[start2:end2]
    return  imgarr
 
def test_writetotif(dir):
    #读取二进制数据并转换成int16类型的数组
    #dir = r"J:\G\workk\2023\zy\LGC\MODIS\MOD09GA_A2004107.sur_refl.int"
    f = open(dir, 'rb')
    fint = np.fromfile(f,dtype = np.int16)
 
    #数据提取
    bands, rows, col =6,  2720,3200
    #imgarr = read_as_bil(fint, bands, rows, col)
    imgarr = read_as_bsq(fint, bands, rows, col)
    #imgarr = read_as_bip(fint, bands, rows, col)
    #将提取的数组存储为tif格式图像.
    #注意这里未设置地理坐标和仿射变换信息，所以不能用ENVI等软件打开。
    savedir=dir.split(".")[0]+"_sur_refl_new.tif"
    #savedir = r"E:\datasets\stfdatasets\Coleambally_Irrigation_Area\CIA\MODIS\2001_281_08oct\MOD09GA_A2001281.sur_refl.tif"
    datatype = gdal.GDT_UInt16
    bands, high, width = imgarr.shape
    driver = gdal.GetDriverByName("GTiff")
    datas = driver.Create(savedir, col, rows, bands, datatype)
 
    #设置地理坐标和仿射变换信息
    """
    image_geotrans=[378000,25,0,6170000,0,25]
    image_projection="AMG55"
    datas.SetGeoTransform(image_geotrans)
    datas.SetProjection(image_projection)"""
    for i in range(bands):
        datas.GetRasterBand(i + 1).WriteArray(imgarr[i])
    del datas
    
 
if __name__ =="__main__":
    in_dir = r'J:\G\workk\2023\zy\LGC\MODIS' # input dir
    
    file_list = os.listdir(in_dir)
    for file in file_list:
        if file.endswith('.int'):
            test_writetotif(in_dir+"\\"+file)
    print("save succfully")

以下是CIA_bil格式转tif


"""
Created on Thu Sep 28 14:47:13 2023
@author: Administrator
"""
 
# -*- coding: utf-8 -*-
import numpy as np
from osgeo import gdal
import os
from pyproj import Proj,transform
from osgeo import osr
 
 
# 依据BSQ存储规则，按照存储完单波段整幅图像后再存储下一波段的存储方法进行提取并存入数组。
def read_as_bsq(dataarr,bands, rows, col):
    imgarr = np.zeros((bands,rows,col))
    for b in range(bands):              #取出每个波段
        start = b * rows * col
        end = start + rows * col
        arrtem = dataarr[start:end]
        for r in range(rows):           #一维数组按行取出，存入对应三维数组。
            start2 = r * col
            end2 = start2 + col
            imgarr[b,r,:] = arrtem[start2:end2]
    return  imgarr
 
def ReadBilFile(bil):
 
    gdal.GetDriverByName('EHdr').Register()
    img = gdal.Open(bil)
    band = img.GetRasterBand(1)
    data = band.ReadAsArray()
    return data
# 依据BIP存储规则，按照一个像素所有波段进行存储完，再存储下一个像素所有波段的存储方法进行提取并存入数组。
def read_as_bip(dataarr,bands, rows, col):
    imgarr = np.zeros((bands,rows,col))
    for r in range(rows):               #按行列遍历每个像元
        for c in range(col):
            if r == 0 :
                pix = c
            else:
                pix = r * col + c
            start = pix * bands
            end = start + bands
            arrtem = dataarr[start:end] #从一维数组中取出每个像元的全波段元素（6个）
            for b in range(bands):
                imgarr[b,r,c] = arrtem[b] # 赋值给对应数组
    return  imgarr
# 依据BIL存储规则，按照存储完一行的所有波段再存储下一行，进行提取并存入数组。
def read_as_bil(dataarr,bands, rows, col):
    imgarr = np.zeros((bands,rows,col))
    for r in range(rows): #取出一行的所有波段
        start = r * col * bands
        end = start + col * bands
        arrtem = dataarr[start:end]
        for b in range(bands): #取出每个波段
            start2 = b * col
            end2 = start2 + col
            imgarr[b,r,:] = arrtem[start2:end2]  #存入数组对应位置
    return  imgarr
def test_writetotif(dir):
    #读取二进制数据并转换成int16类型的数组
    #dir = r"J:\G\workk\2023\zy\LGC\MODIS\MOD09GA_A2004107.sur_refl.int"
    f = open(dir, 'rb')
    fint = np.fromfile(f,dtype = np.int16)
 
    #数据提取
    bands, rows, col =6, 2040, 1720
    imgarr = read_as_bil(fint, bands, rows, col)
    #imgarr = read_as_bsq(fint, bands, rows, col)
    #imgarr = read_as_bip(fint, bands, rows, col)
    #将提取的数组存储为tif格式图像.
    #注意这里未设置地理坐标和仿射变换信息，所以不能用ENVI等软件打开。
    savedir=dir.split(".")[0]+"_sur_refl.tif"
    #savedir = r"E:\datasets\stfdatasets\Coleambally_Irrigation_Area\CIA\MODIS\2001_281_08oct\MOD09GA_A2001281.sur_refl.tif"
    datatype = gdal.GDT_UInt16
    bands, high, width = imgarr.shape
    driver = gdal.GetDriverByName("GTiff")
    datas = driver.Create(savedir, col, rows, bands, datatype)
 
    #设置地理坐标和仿射变换信息
    image_geotrans=[378000,25,0,6170000,0,25]
    #image_projection="20255"
    #image_projection="agd66"
    #agd66 = Proj(init='EPSG:20255')
    #p = Proj(init="EPSG:32650")
    datas.SetGeoTransform(image_geotrans)
    #datas.SetProjection(agd66)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(20255)
    datas.SetProjection(outRasterSRS.ExportToWkt())
    for i in range(bands):
        datas.GetRasterBand(i + 1).WriteArray(imgarr[i])
    del datas
    
 
if __name__ =="__main__":
    in_dir = r'J:\G\workk\2023\zy\CIA\Landsat-init' # input dir
    
    file_list = os.listdir(in_dir)
    for file in file_list:
        if file.endswith('.bil'):
            test_writetotif(in_dir+"\\"+file)
    print("save succfully")

以下是CIA int格式转tif


"""
Created on Thu Sep 28 14:47:13 2023
@author: Administrator
"""
 
# -*- coding: utf-8 -*-
import numpy as np
from osgeo import gdal
import os
from osgeo import osr
 
# 依据BSQ存储规则，按照存储完单波段整幅图像后再存储下一波段的存储方法进行提取并存入数组。
def read_as_bsq(dataarr,bands, rows, col):
    imgarr = np.zeros((bands,rows,col))
    for b in range(bands):              #取出每个波段
        start = b * rows * col
        end = start + rows * col
        arrtem = dataarr[start:end]
        for r in range(rows):           #一维数组按行取出，存入对应三维数组。
            start2 = r * col
            end2 = start2 + col
            imgarr[b,r,:] = arrtem[start2:end2]
    return  imgarr
 
def test_writetotif(dir):
    #读取二进制数据并转换成int16类型的数组
    #dir = r"J:\G\workk\2023\zy\LGC\MODIS\MOD09GA_A2004107.sur_refl.int"
    f = open(dir, 'rb')
    fint = np.fromfile(f,dtype = np.int16)
 
    #数据提取
    bands, rows, col =6, 2040, 1720
    #imgarr = read_as_bil(fint, bands, rows, col)
    imgarr = read_as_bsq(fint, bands, rows, col)
    #imgarr = read_as_bip(fint, bands, rows, col)
    #将提取的数组存储为tif格式图像.
    #注意这里未设置地理坐标和仿射变换信息，所以不能用ENVI等软件打开。
    savedir=dir.split(".")[0]+"_sur_refl.tif"
    #savedir = r"E:\datasets\stfdatasets\Coleambally_Irrigation_Area\CIA\MODIS\2001_281_08oct\MOD09GA_A2001281.sur_refl.tif"
    datatype = gdal.GDT_UInt16
    bands, high, width = imgarr.shape
    driver = gdal.GetDriverByName("GTiff")
    datas = driver.Create(savedir, col, rows, bands, datatype)
 
    #设置地理坐标和仿射变换信息
    image_geotrans=[378000,25,0,6170000,0,25]
    image_projection="AMG55"
    datas.SetGeoTransform(image_geotrans)
    
    #datas.SetProjection(agd66)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(20255)
    datas.SetProjection(outRasterSRS.ExportToWkt())
    for i in range(bands):
        datas.GetRasterBand(i + 1).WriteArray(imgarr[i])
    del datas
    
 
if __name__ =="__main__":
    in_dir = r'J:\G\workk\2023\zy\CIA\MODIS-init' # input dir
    
    file_list = os.listdir(in_dir)
    for file in file_list:
        if file.endswith('.int'):
            test_writetotif(in_dir+"\\"+file)
    print("save succfully")

以下是多波段tif数据选择rgb波段然后显示的程序


"""
Created on Sun Oct  8 15:09:42 2023
@author: Administrator
"""
import numpy as np
import tifffile as tf           #tifffile是tiff文件的读取库
from PIL import Image
import cv2 as  cv
 
from osgeo import gdal
import matplotlib.pyplot as plt
#以下为三种拉伸方式，如果不拉伸，图像太黑，拉伸完显示的图像更好看
def optimized_linear(arr):
    a, b = np.percentile(arr, (2.5, 99))
    c = a - 0.1 * (b - a)
    d = b + 0.5 * (b - a)
    arr = (arr - c) / (d - c) * 255
    arr = np.clip(arr, 0, 255)
    return np.uint8(arr)
def percent_linear(arr, percent=2):
    arr_min, arr_max = np.percentile(arr, (percent, 100-percent))
    arr = (arr - arr_min) / (arr_max - arr_min) * 255
    arr = np.clip(arr, 0, 255)
    return np.uint8(arr)
def linear(arr):
    arr_min, arr_max = arr.min(), arr.max()
    arr = (arr - arr_min) / (arr_max - arr_min) * 255
    arr = np.clip(arr, 0, 255)
    return np.uint8(arr)
path=r"J:\G\workk\2023\zy\LGC\MODIS\MOD09GA_A2004235_sur_refl.tif"
data = gdal.Open(path)   # 读取tif文件
num_bands = data.RasterCount     # 获取波段数
print(num_bands)
tmp_img = data.ReadAsArray()      #将数据转为数组
print(tmp_img.shape)
img_rgb = tmp_img.transpose(1, 2, 0)     #由波段、行、列——>行、列、波段
 
img_rgb = np.array(img_rgb, dtype=np.uint16)   #设置数据类型，np.unit8可修改
#img_rgb = np.array(img_rgb) 
r = img_rgb[:, :,3]
g = img_rgb[:, :,2]
b =img_rgb[:, :,1]
 
img_rgb = np.dstack((percent_linear(r), percent_linear(g),percent_linear(b)))    # 波段组合
 
#img_rgb = np.array(img_rgb, dtype=np.uint8) 
 
# plt.imshow(img_rgb)
# plt.show()
# 通过调用plt.axis（“ off”），可以删除编号的轴
plt.figure(dpi=800)
plt.axis("off")
plt.imshow(img_rgb)
plt.show()

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/从前慢现在也慢/article/detail/322843