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()
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
