Sentinel 2 Atmospheric Correction in Google Earth Engine(使用GEE进行哨兵2数据大气校正)
2017-12-28 22:13
7273 查看
From and thanks to:
github/samsammurphy/gee-atmcorr-S2
Import modules
In [3]:
import ee from Py6S import * import datetime import math import os import sys sys.path.append(os.path.join(os.path.dirname(os.getcwd()),'bin')) from atmospheric import Atmospheric ee.Initialize()
time and place
Define the time and place that you are looking for.In [4]:
date = ee.Date('2017-01-01') geom = ee.Geometry.Point(-157.816222, 21.297481) geom = ee.Geometry.Point(12.670733, 41.826685)# ESRIN (ESA Earth Observation Centre)
an image
The following code will grab the first scene that occurs on or after date.In [5]:
# The first Sentinel 2 image S2 = ee.Image( ee.ImageCollection('COPERNICUS/S2') .filterBounds(geom) .filterDate(date,date.advance(3,'month')) .sort('system:time_start') .first() ) # top of atmosphere reflectance toa = S2.divide(10000)
metadata
In [11]:info = S2.getInfo()['properties'] scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start']/1000)# i.e. Python uses seconds, EE uses milliseconds solar_z = info['MEAN_SOLAR_ZENITH_ANGLE']
atmospheric constituents
In [12]:h2o = Atmospheric.water(geom,date).getInfo() o3 = Atmospheric.ozone(geom,date).getInfo() aot = Atmospheric.aerosol(geom,date).getInfo()
target altitude (km)
In [13]:SRTM = ee.Image('CGIAR/SRTM90_V4')# Shuttle Radar Topography mission covers *most* of the Earth alt = SRTM.reduceRegion(reducer = ee.Reducer.mean(),geometry = geom.centroid()).get('elevation').getInfo() km = alt/1000 # i.e. Py6S uses units of kilometers
6S object
The backbone of Py6S is the 6S (i.e. SixS) class. It allows you to define the various input parameters, to run the radiative transfer code and to access the outputs which are required to convert radiance to surface reflectance.In [14]:
# Instantiate s = SixS() # Atmospheric constituents s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o,o3) s.aero_profile = AeroProfile.Continental s.aot550 = aot # Earth-Sun-satellite geometry s.geometry = Geometry.User() s.geometry.view_z = 0 # always NADIR (I think..) s.geometry.solar_z = solar_z # solar zenith angle s.geometry.month = scene_date.month # month and day used for Earth-Sun distance s.geometry.day = scene_date.day # month and day used for Earth-Sun distance s.altitudes.set_sensor_satellite_level() s.altitudes.set_target_custom_altitude(km)
Spectral Response functions
Py6S uses the Wavelength class to handle the wavelength(s) associated with a given channel (a.k.a. waveband). This might be a single scalar value (e.g. a central wavelength) or, if known, possibly the spectral response function of the waveband. The Sentinel2 spectral response functions are provided with Py6S (as well as those of a number of missions). For more details please see the docsor
the (comment-rich) source code
In [15]:
def spectralResponseFunction(bandname): """ Extract spectral response function for given band name """ bandSelect = { 'B1':PredefinedWavelengths.S2A_MSI_01, 'B2':PredefinedWavelengths.S2A_MSI_02, 'B3':PredefinedWavelengths.S2A_MSI_03, 'B4':PredefinedWavelengths.S2A_MSI_04, 'B5':PredefinedWavelengths.S2A_MSI_05, 'B6':PredefinedWavelengths.S2A_MSI_06, 'B7':PredefinedWavelengths.S2A_MSI_07, 'B7':PredefinedWavelengths.S2A_MSI_07, 'B8A':PredefinedWavelengths.S2A_MSI_09, 'B9':PredefinedWavelengths.S2A_MSI_10, 'B10':PredefinedWavelengths.S2A_MSI_11, 'B11':PredefinedWavelengths.S2A_MSI_12, 'B12':PredefinedWavelengths.S2A_MSI_13, } return Wavelength(bandSelect[bandname])
TOA Reflectance to Radiance
Sentinel 2 data is provided as top-of-atmosphere reflectance. Lets convert this to at-sensor radiance for the atmospheric correction.**You can atmospherically corrected directly from TOA reflectance. However, I suggest
radiance for a couple of reasons. Firstly, it is more intuitive. Instead of spherical albedo (which I suspect is more of a mathematical convenience than a physical property) you can use solar irradiance, transmissivity,
path radiance, etc. Secondly, Py6S seems to be more geared towards converting from radiance to SR</sup>
In [16]:
def toa_to_rad(bandname): """ Converts top of atmosphere reflectance to at-sensor radiance """ # solar exoatmospheric spectral irradiance ESUN = info['SOLAR_IRRADIANCE_'+bandname] solar_angle_correction = math.cos(math.radians(solar_z)) # Earth-Sun distance (from day of year) doy = scene_date.timetuple().tm_yday d = 1 - 0.01672 * math.cos(0.9856 * (doy-4))# http://physics.stackexchange.com/questions/177949/earth-sun-distance-on-a-given-day-of-the-year # conversion factor multiplier = ESUN*solar_angle_correction/(math.pi*d**2) # at-sensor radiance rad = toa.select(bandname).multiply(multiplier) return rad
Radiance to Surface Reflectance
Reflected sunlight can be described as follows (wavelength dependence is implied):where L is at-sensor radiance,
is
transmissivity,
is
surface reflectance,
is
direct solar irradiance,
is
diffuse solar irradiance and
is
path radiance. There are five unknowns in this equation, 4 atmospheric terms (
,
,
and
)
and surface reflectance. The 6S radiative transfer code is used to solve for the atmospheric terms, allowing us to solve for surface reflectance.
In [17]:
def surface_reflectance(bandname): """ Calculate surface reflectance from at-sensor radiance given waveband name """ # run 6S for this waveband s.wavelength = spectralResponseFunction(bandname) s.run() # extract 6S outputs Edir = s.outputs.direct_solar_irradiance #direct solar irradiance Edif = s.outputs.diffuse_solar_irradiance #diffuse solar irradiance Lp = s.outputs.atmospheric_intrinsic_radiance #path radiance absorb = s.outputs.trans['global_gas'].upward #absorption transmissivity scatter = s.outputs.trans['total_scattering'].upward #scattering transmissivity tau2 = absorb*scatter #total transmissivity # radiance to surface reflectance rad = toa_to_rad(bandname) ref = rad.subtract(Lp).multiply(math.pi).divide(tau2*(Edir+Edif)) return ref
Atmospheric Correction
In [18]:# original # b = surface_reflectance('B2') # g = surface_reflectance('B3') # r = surface_reflectance('B4') # rgb = r.addBands(g).addBands(b) # all wavebands output = S2.select('QA60') for band in ['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12']: print(band) output = output.addBands(surface_reflectance(band))
B1 B2 B3 B4 B5 B6 B7 B8
--------------------------------------------------------------------------- KeyError Traceback (most recent call last) <ipython-input-18-0eaa6deade39> in <module>() 9 for band in ['B1','B2','B3','B4','B5','B6','B7','B8','B8A','B9','B10','B11','B12']: 10 print(band) ---> 11 output = output.addBands(surface_reflectance(band)) 12 <ipython-input-17-31b567fffbcd> in surface_reflectance(bandname) 5 6 # run 6S for this waveband ----> 7 s.wavelength = spectralResponseFunction(bandname) 8 s.run() 9 <ipython-input-15-46be73bfbf09> in spectralResponseFunction(bandname) 19 } 20 ---> 21 return Wavelength(bandSelect[bandname]) KeyError: 'B8'
Display results
In [ ]:from IPython.display import display, Image region = geom.buffer(5000).bounds().getInfo()['coordinates'] channels = ['B4','B3','B2'] original = Image(url=toa.select(channels).getThumbUrl({ 'region':region, 'min':0, 'max':0.25 })) corrected = Image(url=ref.select(channels).getThumbUrl({ 'region':region, 'min':0, 'max':0.25 })) display(original, corrected)
Export to Asset
In [ ]:# # set some properties for export # dateString = scene_date.strftime("%Y-%m-%d") # ref = ref.set({'satellite':'Sentinel 2', # 'fileID':info['system:index'], # 'date':dateString, # 'aerosol_optical_thickness':aot, # 'water_vapour':h2o, # 'ozone':o3})
In [ ]:
# define YOUR assetID # in my case it was something like this.. # assetID = 'users/samsammurphy/shared/sentinel2/6S/ESRIN_'+dateString
In [ ]:
# # export # export = ee.batch.Export.image.toAsset(\ # image=ref, # description='sentinel2_atmcorr_export', # assetId = assetID, # region = region, # scale = 30) # # uncomment to run the export # export.start()
相关文章推荐
- 使用Google Earth Engine(一):提取下载MODIS、Landsat点数据
- 使用GDAL对HDF数据进行geoloc校正
- 使用GDAL工具对FY3系列卫星数据进行校正
- 使用Google App Engine进行软件的开发和部署发布
- Google Earth Engine(GEE)学习笔记 零
- 关于Google Earth Engine(GEE)学习笔记搬家说明
- 使用GDAL工具对FY3系列卫星数据进行校正
- 使用GDAL对HDF数据进行校正
- sentinel-2数据下载 大气校正 转ENVI格式
- [转]使用Spatial Adjustment工具进行空间数据校正
- 使用GDAL对HDF数据进行校正
- 使用GDAL对HDF数据进行geoloc校正
- oracle 使用instr()函数对in查询出的记录按照in中的数据进行排序
- Google Earth Engine(GEE)学习笔记 一
- Google Earth Engine(GEE)学习笔记 三
- 使用GDAL对静止卫星圆盘数据进行校正(以FY2为例子)
- 使用 Google App Engine 实现基于云计算的小型 Java 数据服务应用
- 在.net中使用google appengine datastore数据存储(a Soap wrapper)