发布时间:2024-06-24 13:01
目录
1. 实验内容
反权重插值算法简介:
我实现了什么
2. 编程
代码部分
部分运行结果显示(与Arcgis对比):
编辑
使用Arcgiss打开PM2.5 tiff文件(原本图像时黑白变化——》现经过了一些颜色处理)
这是Arcgis进行反距离权重插值的参数设置(为确保与代码结果能够比较,这里将输出像元大小更改为0.001,幂为2默认, 选取的点数更改为11):
这是Arcgis的的结果(颜色也经过了一些处理):
这是ENVI_IDL代码的结果:
使用反距离权重插值法选取距离最近的n个点进行插值
另外这里面用到了新函数sort(),但是应该不难,可以看ENVI的help,也可以自己敲代码熟悉一下这个函数。
pro week_six_test
; 使用反距离权重插值法选取距离最近的n个点进行插值
; 路径
in_dir = 'D:/IDL_program/experiment_data/chapter_4/air_quality_data.csv'
out_dir = 'D:/IDL_program/experiment_data/chapter_4/'
; 获取数据(csv里面第0、1列是经纬度信息、后面几列是诸如pm2.5、pm10等等数据)
data = read_csv(in_dir, header=par_name)
; 获取列索引的相关信息
data_index_name = strsplit(par_name, /extract)
data_index_count = n_elements(data_index_name)
; 提取lon、lat数据
lon = data.(0)
lat = data.(1)
; 假定输出的分辨率
resolution = 0.001 ; 其实分辨率是多少都无所谓的,取决于你需要多大分辨率,因为值都是预测出来的
; 寻找最大最小的经纬度
lon_min = min(lon)
lon_max = max(lon)
lat_min = min(lat)
lat_max = max(lat)
; 计算构建数组的行列数
column = ceil((lon_max - lon_min) / resolution)
row = ceil((lat_max - lat_min) / resolution)
; 计算每个像元的行列号
column_pos = floor((lon - lon_min) / resolution)
row_pos = floor((lat_max - lat) / resolution)
; 对data的后面几列数据(pm2.5\NO2\SO2等等等)进行空间插值
for data_index_i = 2, data_index_count - 1 do begin
; 时间开始计入
start_time = systime(1)
; 输出路径
out_path = out_dir + 'air_qulity_' + data_index_name[data_index_i] + '.tiff'
; 创建存储数据的数组
data_box = fltarr(column, row)
; 放入已知数据
data_box[column_pos, row_pos] = data.(data_index_i)
; 已知数据(样本值)的个数
data_count = n_elements(data.(data_index_i))
; 创建存储预测数据(输出数据)的数组
data_box_out = fltarr(column, row)
; 使用反权重插值公式进行预测数据
for column_i = 0, column - 1 do begin
for row_i = 0, row - 1 do begin
; 首先判断目前值是为已知值(即不是0.0值),如果是那么原样输出即可
if data_box[column_i, row_i] eq 0.0 then begin
distance_sum = 0.0
value_sum = 0.0
; 目前该像元点的经纬度坐标(计算距离需要使用这些)
temp_lon = lon_min + column_i * resolution
temp_lat = lat_max - row_i * resolution
; 对已知样本值进行循环——》计算distance(也就是公式里面的Di的求和)
; 下面的方法是取所有点进行公式的运算
; for data_i = 0, data_count - 1 do begin
; Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
; distance_sum += 1.0 / (Di ^ 2.0) ; 这里公式里面的p我们一般取2,这里也是
; endfor
; for data_i = 0, data_count - 1 do begin
; Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
; value_sum += (1.0 / (Di ^ 2.0)) * data.(data_index_i)[data_i] / distance_sum
; endfor
; 下面的方法是可以选取最近的n个点进行计算
; 由于我们可能会说有这么一个需求:取最近的n个点进行计算(这里会使用sort()函数)
; 由于我们获取最近的n个点,那么就需要所有点的距离,然后排序才知道
; 创建存储每个样本值所在像元与所求像元距离的数组
distance_array = fltarr(data_count)
; 进入循环求得每一个样本值所在像元与所求像元地距离并存储
for data_i = 0, data_count - 1 do begin ; 还有就是下面的平方什么的最好都加上.0表示小数,否则有可能最后会以整数进行运算——》导致结果的小数部分被舍去
Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
distance_array[data_i] = Di
endfor
; 然后对distance_array数组进行排序——》使用sort()函数,传入数组distance_array
near_distance_index = sort(distance_array)
; 第一,sort()函数不会对传入的数组进行更改,第二,返回的也是一个数组,这个数组是传入数组中的元素按从小到大(元素值)排列的索引(实在不理解可以自己看help,也可以自己写几行代码看看是怎么回事(我就是这么干的))
; 这里最近距离的点的个数n的设置
near_n = 11 ; 假如我们设置为11(当然你可以修改设置为其它值)
; 进行取near_n = 11的distance_sum的计算
for near_i = 0, near_n - 1 do begin
data_i = near_distance_index[near_i]
Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
distance_sum += 1.0 / (Di ^ 2.0) ; 公式里面的p我这里选取默认2,(Arcgis里面默认是2)
endfor
; 进行最终值(value_sum)的运算
for near_i = 0, near_n - 1 do begin
data_i = near_distance_index[near_i]
Di = sqrt((temp_lon - lon[data_i]) ^ 2.0 + (temp_lat - lat[data_i]) ^ 2.0)
value_sum += (1.0 / (Di ^ 2.0)) * data.(data_index_i)[data_i] / distance_sum
endfor
; 输出预测的数据
data_box_out[column_i, row_i] = value_sum
endif else begin
data_box_out[column_i, row_i] = data_box[column_i, row_i]
endelse
endfor
endfor
; 地理信息geotiff填写
geo_info={$
MODELPIXELSCALETAG:[resolution,resolution,0.0],$
MODELTIEPOINTTAG:[0.0,0.0,0.0,lon_min,lat_max,0.0],$
GTMODELTYPEGEOKEY:2,$
GTRASTERTYPEGEOKEY:1,$
GEOGRAPHICTYPEGEOKEY:4326,$
GEOGCITATIONGEOKEY:'GCS_WGS_1984',$
GEOGANGULARUNITSGEOKEY:9102,$
GEOGSEMIMAJORAXISGEOKEY:6378137.0,$
GEOGINVFLATTENINGGEOKEY:298.25722}
; 将data中该列的数据输出为geotiff文件
write_tiff, out_path, data_box_out, geotiff=geo_info, /float
stop_time = systime(1)
print, data_index_name[data_index_i] + '>>> ' + strcompress(string(stop_time - start_time)) + 's'
endfor
print, '>>> 程序运行完毕'
end
可以发现,两者几乎没有区别,这说明我们的代码是有效的@!
———————————————————————————————————————————
———————————————————————————————————————————
我是炒茄子,不用谢
MIT 曝光 M1 不可修复的漏洞,网友:时间挺巧,是时候换 M2 了!
(九)PyTorch深度学习:卷积神经网络( GoogleNet网络架构中的 inception module 模块为本次卷积神经网络架构)
Uniapp打开外部网页、跳转指定链接 (打开app以外的web页面)
笔记【Two-Stream Convolutional Networks for Action Recognition in Videos】
Failed to Create Hive Metastore Database Tables.