“””
暗像元校正
方法：计算波段中的最小值进行暗像元校正
过程：1.影像DN值转为辐射能量值（L = DN*gain +bias）
2.计算相对反射率（ρ=π×D2×L/（ESUNI×COS(SZ)））
3.选择波段最小值进行波段值运算（dataset.data =bands-bandsMin）
4.新建栅格数据写入数据（driver = gdal.Create(self, *args, **kwargs)）
@time：2018/4/2 9:40
@author：wujd
“””
from osgeo import gdal
import numpy as np
import os
打开数据文件
os.chdir(r”F:\Wujd\0330darkPixel\3285750”)
dataset = gdal.Open(“HJ1A-CCD1-451-63-20171202-L20003285750-1.TIF”)
获取栅格数据的行列数，波段数，放射矩阵，投影信息
im_width = dataset.RasterXSize
im_height = dataset.RasterYSize
im_bands = dataset.RasterCount
im_geotrans = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
im_data = dataset.ReadAsArray(0,0,im_width,im_height)
type(im_data),im_data.shape
print(im_width, im_height,im_bands,im_geotrans,im_proj,im_data)
3.roi裁切后：遍历部分区域像元，获取最小值
os.chdir(r”F:\Wujd\0330darkPixel”)
datasetRoi = gdal.Open(“roi_min.tif”)
bandRoi = datasetRoi.GetRasterBand(1)
minValue = bandRoi.ComputeRasterMinMax()[0]
print(minValue)
1.计算DN转为辐射能量值——HJ1A-CCD1 b1：gain 1.4609 Bias 7.325 图像的增益与偏置
gain =1.4609
bias = 7.3250
band = dataset.GetRasterBand(1)
data = band.ReadAsArray(0,0,im_width,im_height)
outband_L = (data-minValue)*gain +bias
print(“辐射能量值\n”,outband_L)
print(band.ComputeRasterMinMax())
2.计算相对反射率 ρ=π×D2×L/（ESUNI×COS(SZ)）
days = 237 #拍摄卫片的日期在那一年的天数
sunEle = 21.945 #太阳高度角
ESUNI = 1914.324 #大气顶层太阳辐照度（ESUNI）
D = 1-0.01674*np.sin(np.pi/180*(2*np.pi*(days-93.5)/360)) #日地天文单位距离D
sunZen = 90.000-sunEle #太阳天顶角
band_ref = np.pi*D*outband_L/(ESUNI*np.cos((np.pi/180)*sunZen)) #相对反射率
print(“相对反射率\n”,band_ref)
3.选则波段最小值
arr = band.ReadAsArray(0,0,im_width,im_height).astype(np.float16)
d = np.where(arr >0 )
print(d)
“”“
遍历所有像元
i ,j= 0,0
print(im_height)
for i in range(0,im_height):
for j in range(0,im_width):
if band1[i][j]==0:
band1[i][j]=None
print(band1[i])
i += 1
j += 1
print(band1)
“”“
4.栅格创建写入数据
driver = gdal.GetDriverByName(“GTiff”)
newdataset = driver.Create(“F://Wujd//0330darkPixel//newfile111.tif”,im_width,im_height,im_bands,6,options=[“INTERLEAVE=PIXEL”])
newdataset.SetGeoTransform(im_geotrans)
newdataset.SetProjection(im_proj)
newdataset.GetRasterBand(1).WriteArray(band_ref*10)