FY-4A天顶角计算-2km

from osgeo import gdal
import numpy as np
import math
from time import *
begin_time = time()

filePath = r'F:\z_test\2km_raw.tif'
dataset = gdal.Open(filePath)
adfGeoTransform = dataset.GetGeoTransform()
im_proj = dataset.GetProjection()
nXSize = dataset.RasterXSize #列数
nYSize = dataset.RasterYSize #行数
im_bands = dataset.RasterCount #波段数
band1 = dataset.GetRasterBand(1)#纬度
data_1 = band1.ReadAsArray(0, 0, nXSize, nYSize)

band2 = dataset.GetRasterBand(2)#经度
data_2 = band2.ReadAsArray(0, 0, nXSize, nYSize)

zenith = np.zeros((nYSize, nXSize))
print(nXSize, nYSize)
for i in range(nYSize):
    for j in range(nXSize):
        if data_1[i][j]<360:
            px_hudu = data_1[i][j] * math.pi / 180  # 纬度
            py_hudu = data_2[i][j] * math.pi / 180  # 经度
            sate_hudu = 104.7 * math.pi / 180
            r = 42200
            R = 6370

            gama = math.acos(math.cos(px_hudu) * math.cos(sate_hudu - py_hudu))
            d = r * math.pow(1 + R * R / (r * r) - 2 * (R / r) * math.cos(gama), 1 / 2)

            zenith[i][j] = (math.asin(r * math.sin(gama) / d)) * 180 / math.pi


driverrr = gdal.GetDriverByName("GTiff")
datasetnew = driverrr.Create(r'F:\z_test\FY_SZA.tif', nXSize, nYSize, 1, gdal.GDT_Float32)
datasetnew.SetGeoTransform(adfGeoTransform)
datasetnew.SetProjection(im_proj)
datasetnew.GetRasterBand(1).WriteArray(zenith)


end_time = time()
run_time = end_time-begin_time
print ('该循环程序运行时间：',run_time)