Matlab把geotif重采样到其他分辨率，保持数据总量变化小（R语言无法做到）

需求：把全球0.083栅格数据重采样到0.5，当然可以把栅格首先转成点再用arcgis插值到0.5，但是效率很低，且不够精确，用matlab的imresize可以快速实现该功能

input_file='F:\Dataset\EARTHSTAT\FertilizerCropSpecific_Geotiff\';
crops={'wheat' 'rice' 'maize' 'soybean' 'barley' 'sorghum' 'millet'  'rapeseed' 'groundnut' 'sunflower'...
'sugarcane' 'potato' 'cassava'  'oilpalm' 'rye' 'sugarbeet'};
Nfer_all=[];
for cc=1:16
CropNH3=imread([input_file 'Fertilizer_' crops{cc} '\' crops{cc} '_NitrogenApplication_Rate.tif']);
CropNH3(isnan(CropNH3)|CropNH3<0)=0;
CropNH3=imresize(CropNH3,[360 720],'bilinear');
CropNH3(isnan(CropNH3)|CropNH3<0)=0;
geotiffwrite(['F:/program/R/Farmsize/data/' crops{cc}  '_NitrogenApplication_Rate.tif'],(CropNH3),R);
end

 

此外，也可用R语言快速实现分辨率的改变，有两种方法：

（1）Resample

rs_NOxdep=raster('K:/Dataset/GC/Dep/yr/TotalNOx_2010.tif')

s <- raster(nrow=2160, ncol=4320)
extent(s)=extent(c(-180,180,-90,90))
rs_NOxdep <- resample(rs_NOxdep,s, method='ngb')

cellStats(rs_NOxdep,"sum")

writeRaster(rs_NOxdep,paste('data/NOydep_',2010,'.tif',sep=""),
            options=c('TFW=YES'), overwrite=TRUE)

（2）aggregate

IASINH3data=stack(paste('data/raster/Attainable_Yield_',crops,'2000.tif',sep=''))*100
IASINH3data <- aggregate(IASINH3data, fact = 0.5/res(IASINH3data)) # aggregate output

注意：这两种转换方式，都不是特别精确，对数据的改动很多，不如matlab的imresize精确，会改变全球排放的budget。R语言的优势是统计，盲目对数据重采样操作会带来较大误差。类似例子也体现在R语言Zonal.status按照polygon求取平均、总量，与ArcGIS计算结果偏差较大