创建马歇尔火灾的GOES-16 ABI气溶胶检测可视化
9287 为本文评分:
3.9
创建马歇尔火灾的GOES-16 ABI气溶胶检测可视化
大卫·斯塔巴克 2022年5月2日，星期一
2021年秋季的科罗拉多州极度干旱，12月30日，当二级飓风级风力袭击该地区时，一场城市野火爆发，导致了科罗拉多州历史上最大的财产损失。马歇尔火灾摧毁了苏必利尔市、路易斯维尔市以及博尔德县未建制地区的一千多座建筑。该地区的居民往往只有几分钟的时间撤离，许多人只能穿着身上的衣服逃离。
除了损失的一千多座建筑外，许多其他房屋和企业也遭受了严重的烟雾损害。虽然极端大风和长期干旱为火灾失控创造了理想条件，但居民们现在面临一个现实：这种情况可能并非罕见的风暴叠加事件。在马歇尔火灾之后，科罗拉多州前岭地区（该地区人口众多，靠近火灾顶级生态系统）正在接受地方性野火和全年皆是的火灾季节这一真实可能性。
利用影像研究野火
地理空间影像可用于研究火灾前的条件、监测火灾期间的情况以及评估灾后损失。来自高级基线成像仪（ABI）的数据（ABI是安装在L3Harris建造的GOES-16平台上的传感器），使我们能够以高时间、空间和光谱分辨率研究所有类型的天气现象。某些波段会突出烟雾，另一些会突出不同大小的气溶胶，还有其他波段则突出火灾本身。能够访问、分析和可视化ABI数据的人越多，我们就能从这类事件中学到更多东西，了解它们如何影响人们的生活、健康、航空安全等。
ABI数据集有直达NOAA以及威斯康星大学等机构的设施用于预报业务，但也通过AWS S3存储桶公开提供。这篇博文提供了代码示例，并展示了任何人如何使用IDL®编程语言从AWS S3中拉取该日期的GOES-16 ABI数据，创建显示烟雾动态和范围的RGB图像，并创建气溶胶检测数据产品。
访问ABI数据的逐步指南
访问ABI数据的第一步是使用IDL中的IDLnetURL对象从AWS S3拉取数据。我使用AWS Amazon S3 API指南来了解如何填写IDLnetURL参数 - Amazon S3 - Amazon Simple Storage Service。根据这个指南，我编写了以下例程来查询Goes-16的S3数据。
function dj_goes_16_s3_query, prefix=prefix, print=print
compile_opt idl2
s3UrlObject = obj_new('idlneturl')
bucket = 'noaa-goes16'
region = 'us-east-1'
s3UrlObject.setProperty, url_host = bucket + '.s3.' + region + '.amazonaws.com'
urlQuery = 'list-type=2&delimiter=/'
if keyword_set(prefix) then begin
urlQuery = urlQuery + '&prefix='+prefix
endif
keyList=list()
dateList = list()
sizeList = list()
gotAll = 0
startAfter = ''
initUrlQuery = urlQuery
while (gotAll eq 0) do begin
urlQuery = initUrlQuery
if startAfter ne '' then begin
urlQuery = urlQuery + '&start-after=' + startAfter
startAfter = ''
endif
s3urlObject.setproperty, url_query= urlQuery
result = s3urlObject.get()
parsedData = xml_parse(result)
if ((parsedData['ListBucketResult']).haskey('Contents')) then begin
contents = ((parsedData['ListBucketResult'])['Contents'])
for awsKeyIndex=0L, n_elements(contents)-1 do begin
keyList.add, (contents[awsKeyIndex])['Key']
dateList.add, (contents[awsKeyIndex])['LastModified']
sizeList.add, (contents[awsKeyIndex])['Size']
endfor
if keyword_set(print) then begin
outArray = strarr(3,n_elements(keyList))
outArray[0,*] = transpose(keyList.toArray())
outArray[1,*] = transpose(dateList.toArray())
outArray[2,*] = transpose(sizeList.toArray())
print, '---------------------------------'
print, 'Key, Date, Size
print, '---------------------------------'
print, outArray
endif
endif
prefixList= list()
if ((parsedData['ListBucketResult']).haskey('CommonPrefixes')) then begin
commonPrefixes = ((parsedData['ListBucketResult'])['CommonPrefixes'])
if n_elements(commonPrefixes) eq 1 then begin
prefixList.add, commonPrefixes['Prefix']
endif else begin
for prefixKeyIndex=0L, n_elements(commonPrefixes)-1 do begin
prefixList.add, (commonPrefixes[prefixKeyIndex])['Prefix']
endfor
endelse
if keyword_set(print) then begin
print, '---------------------------------'
print, 'Prefixes'
print, '---------------------------------'
print, transpose(prefixList.toArray());objListStructor['prefixlist']
endif
endif
if ((parsedData['ListBucketResult'])['IsTruncated'] eq 'true') then begin
startAfter = keyList[-1]
endif else begin
gotAll = 1
endelse
endwhile
returnHash = hash()
returnhash['keylist'] = keyList
returnhash['datelist'] = dateList
returnhash['sizelist'] = sizeList
returnhash['prefixlist'] = prefixList
obj_destroy, s3UrlObject
return, returnhash
end
在此例程中，我使用 "noaa-goes16" 作为存储桶，"us-east-1" 作为区域。
s3UrlObject = obj_new('idlneturl')
bucket = 'noaa-goes16'
region = 'us-east-1'
s3UrlObject.setProperty, url_host = bucket + '.s3.' + region + '.amazonaws.com'
我使用 "list-type=2" 和 "delimiter=/" 作为 "url_query" 属性的参数。运行查询将返回存储桶顶层可用的键（文件）和前缀（子目录）列表。如果没有分隔符，响应会列出存储桶中的所有文件，这将创建一个非常长的列表。
s3urlObject.setproperty, url_query='list-type=2&delimiter=/'
result = s3urlObject.get()
响应将是一个XML对象，可以使用xml_parse例程进行解析。之后，可以提取键和前缀。
parsedData = xml_parse(result)
在没有指定任何前缀的情况下运行 "dj_goes_16_s3_query"，将产生以下结果：
IDL> topLevelQuery = dj_goes_16_s3_query(/print)
---------------------------------
Key, Date, Size
---------------------------------
Beginners_Guide_to_GOES-R_Series_Data.pdf 2021-04-02T17:12:02.000Z 5887019
Version1.1_Beginners_Guide_to_GOES-R_Series_Data.pdf 2021-03-23T15:49:55.000Z 3578396
index.html 2021-09-27T19:48:14.000Z 32357
---------------------------------
Prefixes
---------------------------------
ABI-L1b-RadC/
ABI-L1b-RadF/
ABI-L1b-RadM/
ABI-L2-ACHAC/
ABI-L2-ACHAF/
ABI-L2-ACHAM/
ABI-L2-ACHTF/
ABI-L2-ACHTM/
ABI-L2-ACMC/
ABI-L2-ACMF/
ABI-L2-ACMM/
…
SUVI-L1b-Fe284/
然后，我使用 "prefix" 参数来查询特定前缀下的键和前缀。例如，我使用以下调用来查看 "ABI-L1b-RadC/" 下可用的前缀：
IDL> productTypeLevelQuery = dj_goes_16_s3_query(/print, prefix='ABI-L1b-RadC/')
---------------------------------
Prefixes
---------------------------------
ABI-L1b-RadC/2000/
ABI-L1b-RadC/2017/
ABI-L1b-RadC/2018/
ABI-L1b-RadC/2019/
ABI-L1b-RadC/2020/
ABI-L1b-RadC/2021/
ABI-L1b-RadC/2022/
我继续对这个过程使用进一步的子目录，并为每个键找到了以下结构：
<Product Type>/<year>/<day>/<hour>/<filename.nc>
我参考了GOES-16: True Color Recipe — Unidata Python Gallery教程来了解如何从GOES-16数据生成RGB图像。根据该教程，应使用二级多通道格式化数据产品（ABI-L3-MCIMP）。因此，我使用以下查询来确定2022年12月30日该类型有哪些可用文件：
IDL> hourlevelQuery = dj_goes_16_s3_query(/print, prefix='ABI-L2-MCMIPC/2021/364/17/')
我打印了结果，发现第一个可用的键是：
ABI-L2-MCMIPC/2021/364/17/OR_ABI-L2-MCMIPC-M6_G16_s20213641701172_e20213641703545_c20213641704048.nc
我编写了以下例程来从存储桶下载文件。
function dj_goes_16_s3_download, key, outputdirectory=outputdirectory
compile_opt idl2
s3UrlObject = obj_new('idlneturl')
bucket = 'noaa-goes16'
region = 'us-east-1'
s3UrlObject.setProperty, url_host = bucket + '.s3.' + region + '.amazonaws.com'
fileBaseSplitString = strsplit(key, '/', /EXTRACT)
fileBaseName = fileBaseSplitString[-1]
s3urlObject.setProperty, url_path=Key
if (keyword_set(outputdirectory)) then begin
outputFilename = filepath(filebasename, ROOT_DIR=outputDirectory)
endif else begin
outputFilename = fileBaseName
endelse
;if file already exist, don't re-download it
if (file_test(outputFilename)) then begin
objectFile = outputfilename
endif else begin
objectFile = s3urlObject.get(filename=outputFileName)
endelse
return, objectFile
end
要下载文件，我使用键 "ABI-L2-MCMIPC/2021/364/17/OR_ABI-L2-MCMIPC-M6_G16_s20213641701172_e20213641703545_c20213641704048.nc" 运行了 "dj_goes_16_s3_download"：
IDL>objectFile = dj_goes_16_s3_download('ABI-L2-MCMIPC/2021/364/17/OR_ABI-L2-MCMIPC-M6_G16_s20213641701172_e20213641703545_c20213641704048.nc')
GOES-16数据是NetCDF格式，可以使用ncdf_parse进行解析。
data = ncdf_parse(objectFile, /read_data)
center_lon = data['geospatial_lat_lon_extent', $
'geospatial_lon_nadir', '_DATA']
xscale = data['x', 'scale_factor', '_DATA']
xoffset = data['x', 'add_offset', '_DATA']
x_radians = data['x', '_DATA']*DOUBLE(xscale) + xoffset
yscale = data['y', 'scale_factor', '_DATA']
yoffset = data['y', 'add_offset', '_DATA']
y_radians = data['y', '_DATA']*DOUBLE(yscale) + yoffset
我使用以下命令从数据中获取红色、蓝色和真绿色值。在缩放时，我将最大值设置为2,000，因为我查看了直方图，发现大部分数据都低于该值。
redBand = bytscl(data['CMI_C02', '_DATA'], max=2000)
greenBand = bytscl(data['CMI_C03', '_DATA'], max=2000)
blueBand = bytscl(data['CMI_C01', '_DATA'], max=2000)
trueGreen = bytscl(0.45*redBand + 0.1*greenBand + 0.45*blueBand)
bandDims = size(redBand, /DIMENSIONS)
rgbArray = fltarr(3,bandDims[0],bandDims[1])
rgbArray[0,*,*] = redBand
rgbArray[1,*,*] = trueGreen
rgbArray[2,*,*] = blueBand
然后，我使用image函数来显示数据：
windowDimX = 800
windowDimY = 450
margin = [.2,.2,.1,.1]
i = IMAGE(bytscl(rgbArray), x_radians, y_radians, $
MAP_PROJECTION='GOES-R', $
margin=margin, $
GRID_UNITS='meters', $
CENTER_LONGITUDE=center_lon, $
DIMENSIONS=[windowDimX,windowDimY], $
TITLE='GOES-16 '+ data['time_coverage_start', '_DATA'])
mg = i.MapGrid
mg.label_position = 0
mg.clip = 0
mc = MAPCONTINENTS(/COUNTRIES, COLOR='red')
mc = MAPCONTINENTS(/US, COLOR='red')
然后，我遵循类似的程序来显示气溶胶检测。气溶胶数据是一个包含三个值（0, 1, 255）的字节数组。我添加了一个颜色表，将0设置为黑色，1设置为红色，255设置为白色。
objectFile = dj_goes_16_s3_download('ABI-L2-ADPC/2021/364/17/OR_ABI-L2-ADPC-M6_G16_s20213641701172_e20213641703545_c20213641705442.nc')
data = ncdf_parse(objectFile, /read_data)
center_lon = data['geospatial_lat_lon_extent', $
'geospatial_lon_nadir', '_DATA']
xscale = data['x', 'scale_factor', '_DATA']
xoffset = data['x', 'add_offset', '_DATA']
x_radians = data['x', '_DATA']*DOUBLE(xscale) + xoffset
yscale = data['y', 'scale_factor', '_DATA']
yoffset = data['y', 'add_offset', '_DATA']
y_radians = data['y', '_DATA']*DOUBLE(yscale) + yoffset
aerosol = data['Aerosol', '_DATA']
rgbTable = bytarr(3,255)
rgbTable[*,0] = 255
rgbTable[0,1] = 255
rgbTable[1,1] = 0
rgbTable[2,1] = 0
windowDimX = 800
windowDimY = 450
margin = [.2,.2,.1,.1]
i = IMAGE(aerosol, x_radians, y_radians, $
MAP_PROJECTION='GOES-R', GRID_UNITS='meters', $
margin=margin, $
CENTER_LONGITUDE=center_lon, $
RGB_TABLE=rgbTable, $
DIMENSIONS=[windowDimX,windowDimY], $
TITLE='GOES-16 Aerosol '+ data['time_coverage_start', '_DATA'])
以下程序是一个示例程序，它使用 "dj_goes_16_s3_query" 和 "dj_goes_16_s3_download" 来生成马歇尔火灾期间（UTC时间17:00至22:00）使用GOES-16数据的真彩色图像动画。如果使用气溶胶关键字运行此程序，它将使用气溶胶检测GOES-16数据产品生成动画。
pro dj_example_plot_areosol_goes_aws, aerosol=aerosol
compile_opt idl2
windowDimX = 800
windowDimY = 450
margin = [.2,.2,.1,.1]
topLevelQuery = dj_goes_16_s3_query(/print)
prefixes = (topLevelQuery['prefixlist']).toArray()
multiChannelIndex = where(strmatch(prefixes,'*ABI-L2-MCMIPC*'))
aerosolIndex = where(strmatch(prefixes, '*ABI-L2-ADPC*'))
if keyword_set(aerosol) then begin
productTypeLevelQuery = dj_goes_16_s3_query(prefix=prefixes[aerosolIndex[0]], /print)
endif else begin
productTypeLevelQuery = dj_goes_16_s3_query(prefix=prefixes[multiChannelIndex[0]], /print)
endelse
prefixes = (productTypeLevelQuery['prefixlist']).toArray()
yearIndex = where(strmatch(prefixes,'*2021*'))
yearLevelQuery = dj_goes_16_s3_query(prefix=prefixes[yearIndex[0]], /print)
prefixes = (yearLevelQuery['prefixlist']).toArray()
dayIndex = where(strmatch(prefixes,'*364*'))
dayLevelQuery = dj_goes_16_s3_query(prefix=prefixes[dayIndex[0]], /print)
hoursPrefixes = dayLevelQuery['prefixlist']
;Create video object for projected video
if keyword_set(aerosol) then begin
oVid = IDLffVideoWrite('aerosol_projected_output.mp4',FORMAT='mp4')
vidstream = oVid.AddVideoStream(windowdimx,windowdimy,10)
endif else begin
oVid = IDLffVideoWrite('projected_output.mp4',FORMAT='mp4')
vidstream = oVid.AddVideoStream(windowdimx,windowdimy,10)
endelse
;Store a list of the hours prefixes for
;hoursPrefixes = objListStructor['prefixlist']
for hoursIndex = 17L, 21 do begin
hourLevelQuery = dj_goes_16_s3_query(prefix=hoursPrefixes[hoursIndex], /print)
objectKey = (hourLevelQuery['keylist'])
for objectKeyIndex = 0L, n_elements(objectKey)-1 do begin
objectFile = dj_goes_16_s3_download(objectKey[objectKeyIndex])
data = ncdf_parse(objectFile, /read_data)
print, 'Time start = ', data['time_coverage_start', '_DATA']
; Extract the center longitude and radiance data
center_lon = data['geospatial_lat_lon_extent', $
'geospatial_lon_nadir', '_DATA']
; radiance = data['Rad','_DATA']
; Compute the X/Y pixel locations in radians
xscale = data['x', 'scale_factor', '_DATA']
xoffset = data['x', 'add_offset', '_DATA']
x_radians = data['x', '_DATA']*DOUBLE(xscale) + xoffset
yscale = data['y', 'scale_factor', '_DATA']
yoffset = data['y', 'add_offset', '_DATA']
y_radians = data['y', '_DATA']*DOUBLE(yscale) + yoffset
if (keyword_set(aerosol)) then begin
aerosol = data['Aerosol', '_DATA']
smoke = data['Smoke', '_DATA']
rgbTable = bytarr(3,255)
rgbTable[*,0] = 255
rgbTable[0,1] = 255
rgbTable[1,1] = 0
rgbTable[2,1] = 0
if (obj_valid(i)) then begin
i.setdata, aerosol, x_radians, y_radians
i.title = 'GOES-16 Aerosol '+ data['time_coverage_start', '_DATA']
endif else begin
; Use the GOES-R (GOES-16) map projection
i = IMAGE(aerosol, x_radians, y_radians, $
MAP_PROJECTION='GOES-R', GRID_UNITS='meters', $
margin=margin, $
CENTER_LONGITUDE=center_lon, $
RGB_TABLE=rgbTable, $
DIMENSIONS=[windowDimX,windowDimY], $
TITLE='GOES-16 Aerosol '+ data['time_coverage_start', '_DATA'])
i.limit = [38,-106, 41, -101]
mg = i.MapGrid
mg.label_position = 0
mg.clip = 0
mg.linestyle = 'none'
mg.box_axes = 1
mc = MAPCONTINENTS(/COUNTRIES, COLOR='teal')
mc = MAPCONTINENTS(/US, COLOR='teal')
annotationArrow = arrow([-105,-105.1],[40.5,40.05], /DATA,$
COLOR='blue', $
/CURRENT)
t2 = text(-104.995, 40.51,'Approx. Smoke Origin', color='blue', /DATA, font_style=1)
endelse
time = ovid.put(vidStream, i.copywindow())
print, time
endif else begin
redBand = bytscl(data['CMI_C02', '_DATA'], max=2000)
greenBand = bytscl(data['CMI_C03', '_DATA'], max=2000)
blueBand = bytscl(data['CMI_C01', '_DATA'], max=2000)
trueGreen = bytscl(0.45*redBand + 0.1*greenBand + 0.45*blueBand)
bandDims = size(redBand, /DIMENSIONS)
rgbArray = fltarr(3,bandDims[0],bandDims[1])
rgbArray[0,*,*] = redBand
rgbArray[1,*,*] = trueGreen
rgbArray[2,*,*] = blueBand
;Apply squareroot stetch to improve image clarity.
rgbArray = sqrt(rgbArray/255)*255
; Use the GOES-R (GOES-16) map projection
if obj_valid(i) then begin
;skip frame at 2021-12-30T19:31:17.2Z
if (data['time_coverage_start', '_DATA'] ne '2021-12-30T19:31:17.2Z') then begin
i.setData, rgbArray, x_radians, y_radians
i.title = 'GOES-16 '+ data['time_coverage_start', '_DATA']
endif
endif else begin
i = IMAGE(bytscl(rgbArray), x_radians, y_radians, $
MAP_PROJECTION='GOES-R', $
margin=margin, $
GRID_UNITS='meters', $
CENTER_LONGITUDE=center_lon, $
DIMENSIONS=[windowDimX,windowDimY], $
TITLE='GOES-16 '+ data['time_coverage_start', '_DATA'])
i.limit = [38,-106, 41, -101]
mg = i.MapGrid
mg.label_position = 0
mg.clip = 0
mg.linestyle = 'none'
mg.box_axes = 1
mc = MAPCONTINENTS(/COUNTRIES, COLOR='red')
mc = MAPCONTINENTS(/US, COLOR='red')
annotationArrow = arrow([-105,-105.1],[40.5,40.05], /DATA,$
COLOR='red', $
/CURRENT)
t2 = text(-104.995, 40.51,'Approx. Smoke Origin', color='red', /DATA, font_style=1)
endelse
;skip frame at 2021-12-30T19:31:17.2Z
if (data['time_coverage_start', '_DATA'] ne '2021-12-30T19:31:17.2Z') then begin
time = ovid.put(vidStream, i.copywindow())
print, time
endif
endelse
endfor
endfor
ovid = 0
end
感谢您的关注。请期待更多关于使用GOES、EUMETSAT和其他仪器访问和创建可视化的博客。IDL提供了对复杂数据和代码的灵活访问，并且可以通过IDL虚拟机（一种免费运行IDL代码的方式！）与同事共享。我的团队擅长使用ENVI®和IDL解决图像相关的难题，并连接其他编程语言。如果您有我们可能能够协助的项目，请发送电子邮件至GeospatialInfo@NV5.com。
探索更多
录制的网络研讨会
SAR三部分系列
学习合成孔径雷达（SAR）的应用和优势。
真实案例研究
SAR网络研讨会
使用SAR干涉测量技术进行地表运动监测
灾害管理
录制的网络研讨会
在灾害管理周期中有效使用地理空间数据