240501

数值缺失由数据转换为浮点导致

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
import os # 使用代理
os.environ['HTTP_PROXY'] = 'http://127.0.0.1:10809'
os.environ['HTTPS_PROXY'] = 'http://127.0.0.1:10809'

import ee
ee.Authenticate() # 登录

ee.Initialize(project='ee-dsm-watershed-240315') # 选择项目

import geopandas # 读取本地文件

file = geopandas.read_file(r"C:\Users\Qiao\Documents\RStudio\samp._0430\line_topo.gpkg")

import geemap
line_topo = ee.Feature(geemap.gdf_to_ee(file).geometry(), {'label': 'line_topo'})
Map = geemap.Map(); Map.addLayer(line_topo); Map.centerObject(line_topo, 13); Map

# https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2
inBands = ee.List(['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7'])
outBands = ee.List(['blue','green','red','nir','swir1','swir2'])

# 地形校正代码改编自 https://mygeoblog.com/2018/10/17/terrain-correction-in-gee/

def SCSC(img):
    dem = ee.Image("USGS/SRTMGL1_003"); dem
    scale = 300; degree2radian = 0.01745
    cen_point = img.geometry().centroid(); cen_point
    SZ_rad = (
        # ee.Image.constant(ee.Number(img.get('SOLAR_ZENITH_ANGLE')))
        ee.Image.constant(ee.Number(90.0).subtract(ee.Number(img.get('SUN_ELEVATION'))))
        .multiply(3.14159265359).divide(180)
        .clip(img.geometry().buffer(10000))
    )
    SA_rad = (
        # ee.Image.constant(ee.Number(img.get('SOLAR_AZIMUTH_ANGLE')))
        ee.Image.constant(ee.Number(img.get('SUN_AZIMUTH')))
        .multiply(3.14159265359).divide(180)
        .clip(img.geometry().buffer(10000))
    )
    slp = (
        ee.Terrain.slope(dem)
        .clip(img.geometry().buffer(10000))
    )
    slp_rad = (
        ee.Terrain.slope(dem)
        .multiply(3.14159265359).divide(180)
        .clip(img.geometry().buffer(10000))
    )
    asp_rad = (
        ee.Terrain.aspect(dem)
        .multiply(3.14159265359).divide(180)
        .clip(img.geometry().buffer(10000))
    )
    cosZ = SZ_rad.cos()
    cosS = slp_rad.cos()
    slope_illumination = cosS.expression(
        "cosZ * cosS", 
        {'cosZ': cosZ, 
         'cosS': cosS.select('slope')
        }
    )
    sinZ = SZ_rad.sin()
    sinS = slp_rad.sin()
    cosAziDiff = (SA_rad.subtract(asp_rad)).cos()
    aspect_illumination = sinZ.expression(
        "sinZ * sinS * cosAziDiff", 
        {'sinZ': sinZ,
         'sinS': sinS,
         'cosAziDiff': cosAziDiff
        }
    )
    ic = slope_illumination.add(aspect_illumination)
    img_plus_ic = ee.Image(
        img.addBands(ic.rename('IC'))
        .addBands(cosZ.rename('cosZ'))
        .addBands(cosS.rename('cosS'))
        .addBands(slp.rename('slope'))
    ); img_plus_ic
    props = img_plus_ic.toDictionary(); props
    st = img_plus_ic.get('system:time_start'); st
    mask1 = img_plus_ic.select('nir').gt(-0.1)
    mask2 = (
        img_plus_ic.select('slope').gte(5)
        .And(img_plus_ic.select('IC').gte(0))
        .And(img_plus_ic.select('nir').gt(-0.1))
    )
    img_plus_ic_mask2 = ee.Image(img_plus_ic.updateMask(mask2))
    bandList = ['blue','green','red','nir','swir1','swir2']
    compositeBands = img_plus_ic.bandNames()
    nonCorrectBands = img_plus_ic.select(compositeBands.removeAll(bandList))
    compositeBands;nonCorrectBands
    geom = (
        ee.Geometry(img_plus_ic.get('system:footprint'))
        .bounds().buffer(10000)
    )
    def apply_SCSccorr(band):
        method = 'SCSc'
        out = (
            img_plus_ic_mask2.select('IC', band)
            .reduceRegion(
                reducer = ee.Reducer.linearFit(), 
                geometry = ee.Geometry(img_plus_ic.geometry().buffer(-5000)), 
                scale = 300,
                maxPixels = 1000000000
            )
        )
        if out is None:
            return img_plus_ic_mask2.select(band)
        else:
            out_a = ee.Number(out.get('scale'))
            out_b = ee.Number(out.get('offset'))
            out_c = out_b.divide(out_a)
            SCSc_output = img_plus_ic_mask2.expression(
                "((image * (cosB * cosZ + cvalue)) / (ic + cvalue))", 
                {'image': img_plus_ic_mask2.select(band),
                 'ic': img_plus_ic_mask2.select('IC'),
                 'cosB': img_plus_ic_mask2.select('cosS'),
                 'cosZ': img_plus_ic_mask2.select('cosZ'),
                 'cvalue': out_c
                }
            )
            return SCSc_output
    bandList_SCS = []
    for band in bandList:
        image = apply_SCSccorr(band)
        bandList_SCS.append(image)
    bandList_SCS
    img_SCSccorr = (
        ee.Image(bandList_SCS)
        .addBands(img_plus_ic.select('IC'))
    );img_SCSccorr
    bandList_IC = ee.List([bandList, 'IC']).flatten()
    img_SCSccorr_R = (
        (img_SCSccorr.unmask(img_plus_ic.select(bandList_IC)).select(bandList))
        .addBands(nonCorrectBands).setMulti(props).set('system:time_start',st)
    );img_SCSccorr_R
    # 此处如不加getinfo则为GEE格式，添加后为字符串
    # 注意'system:id'和'id'的区别
    newimg = img_SCSccorr_R.set({'system:id': img_SCSccorr_R.get('L1_LANDSAT_PRODUCT_ID').getInfo()})
    return ee.Image(newimg)

dataset = (
    # https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2
    ee.ImageCollection('LANDSAT/LC09/C02/T1_L2')
    .filterDate('2023-09-20', '2023-11-01')
    .filterBounds(line_topo.geometry().centroid())
    .filter(ee.Filter.lt('CLOUD_COVER_LAND',25))
    .select(inBands,outBands)
);dataset

Map.add_layer(dataset.first(),{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'test');Map
Map.add_layer(SCSC(dataset.first()),{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'test_corr');Map

dataset.first()

# LANDSAT/LC08/C02/T1_L2/LC08_122023_20230314
# LANDSAT/LC09/C02/T1_L2/LC09_123023_20230617
# LANDSAT/LC08/C02/T1_L2/LC08_122023_20230805
# LANDSAT/LC09/C02/T1_L2/LC09_122023_20230930

def apply_scale_factors(image):
  optical_bands = image.select(['blue','green','red','nir','swir1','swir2']).multiply(0.0000275).add(-0.2)
  return image.addBands(optical_bands, None, True)

def NDVI(image): 
    ndvi = image.normalizedDifference(['nir', 'red']).rename('NDVI')
    return image.addBands(ndvi).select(['NDVI'])

img1 = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_122023_20230314').select(inBands,outBands)
img2 = ee.Image('LANDSAT/LC09/C02/T1_L2/LC09_123023_20230617').select(inBands,outBands)
img3 = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_122023_20230805').select(inBands,outBands)
img4 = ee.Image('LANDSAT/LC09/C02/T1_L2/LC09_122023_20230930').select(inBands,outBands)
Map = geemap.Map(); Map.addLayer(line_topo); Map.centerObject(line_topo, 13); Map
Map.add_layer(img1,{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'img1')
Map.add_layer(img2,{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'img2')
Map.add_layer(img3,{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'img3');Map
Map.add_layer(img4,{ 'bands': ['red', 'green', 'blue'], 'min': 0, 'max': 20000},'img4');Map

# 重要！数值缺失由数据转换为浮点导致！！
img1 = ee.Image('LANDSAT/LC08/C02/T1_L2/LC08_122023_20230314').select(inBands,outBands)
img1_SCSC = SCSC(img1)
img1_SCSC_NDVI = NDVI(img1_SCSC)
Map = geemap.Map();Map.centerObject(line_topo, 13)
Map.add_layer(img1_SCSC_NDVI,{ 'bands': 'NDVI', 'min': -1, 'max': 1},'NDVI');Map

geemap.download_ee_image(
    img1_SCSC_NDVI, 
    filename = r"C:\Users\Qiao\Documents\RStudio\samp._0430\img1_SCSC_NDVI.tif", 
    region = geemap.gdf_to_ee(file).geometry()
)