01 要求
获取一张地表温度数据,对其进行分组统计,并制作间隔为3K的温度分布频数和频率直方图。
02 下载地表温度数据
首先需要下载地表温度数据,我想了几种产品,一个是MODIS的MOD11A2产品,一个是Landsat8/9 的C2,还有sentinel3的SLSTR传感器提供的地表温度数据,当然后面查阅发现还有VIIRS、ASTER(Terra卫星)等等,这里不展开说明。
由于MODIS的产品下载即可使用,无需过多的预处理,这里从NASA下载MOD11A2产品吧。
2.1 检索
好吧,下载的居然还是HDF文件,为了不偏题(没时间),这里就不展开说明如何得到tiff文件。
本来是打算从头到尾讲述一下,这么一下工作量大了些。
如下是得到的tiff文件:
03 直方图绘制
3.1 IDL代码
说明都在注释中,如下:
; 该程序用于计算某幅影像中各个分组下的频率/频数以及分布直方图
pro histogram_plot
img_path = 'E:\FeaturesTargets\uniform\LST\LST_2000_03.tiff' ; 输入路径
img = read_tiff(img_path) ; 读取栅格矩阵, float[139列, 184行]
; 设置无效值为nan, 此处-65535.000为无效值
img[where(img eq -65535.000)] = !values.f_nan
img_min = min(img, /nan) ; 取栅格矩阵的最小值, /nan表示忽略nan
img_max = max(img, /nan) ; 取栅格矩阵的最大值, /nan同上
; 获取各个范围的频数分布(实际上对于这种不是特别特殊的图建议使用excel绘制会快很多)
hist = histogram(img, binsize=3.0, /nan, min=img_min, max=img_max) ; 计算各个范围的频数
x_range = findgen(ceil((img_max - img_min) / 3)) * 3 + img_min ; 得到x范围
fig = barplot(x_range, hist, /histogram, xrange=[img_min, img_max], xtickinterval=3.0, $
xtitle='温度范围', ytitle='温度(频数)', title='某地温度频数分布直方图', font_name='Microsoft Yahei', font_size=14)
fig.save, 'D:\Objects\JuniorSpringTerm\GeoMathAnalysis\Chapter2\histogram1.png'
; 上述是频数直方图绘制, 下面是频率直方图绘制, 频率等于频数乘以总数
valid_sum = total(~finite(img, /nan)) ; 或者total(hist)
hist = hist / valid_sum
fig = barplot(x_range, hist, /histogram, xrange=[img_min, img_max], xtickinterval=3.0, $
xtitle='温度范围', ytitle='温度(频数)', title='某地温度频率分布直方图', font_name='Microsoft Yahei', font_size=14)
fig.save, 'D:\Objects\JuniorSpringTerm\GeoMathAnalysis\Chapter2\histogram2.png'
end
其中涉及的两个比较关键的函数为histogram
和barplot
。
对于histogram
,其是计算各个分组下的频数(以一维数组形式)。img
表示传入的栅格矩阵(数组),binsize
表示分组间隔,min
表示进行分组的最小值,max
表示进行分组的最大值,数组中不在最大值和最小值范围之间的数不纳入计算频数(当然你也可以不传入,默认会计算传入数组的最大最小值)。
对于barplot
,这就是进行条形图及其变种图的绘制,这里计算到这里你也可以大可不用代码进行绘制如果没有自动化和大批量的要求的话。x_range
表示传入的X,hist表示传入的Y,/histogram
表示按照直方图形式绘制,xrange
表示x轴绘制的最大最小值范围,xtickinterval
表示X轴各个数值标记或者说绘制的间隔,xtitle
表示X轴标题,ytitle
表示Y轴标题,title
表示图主标题。