Matlab求取geotiff不同排放国家均值、总量

1、求取0.5°栅格面积（R语言）

这个求取全球0.5°栅格像元面积matlab可能也可以计算，但是比较麻烦，R语言有自带的包可以计算，输出成geotiff

library(maptools)
library(raster)
library(marmap)

s <- raster(nrow=360, ncol=720)
extent(s) <-c(-180,180,-90,90)
crs(s) <- CRS("+proj=longlat +datum=WGS84") 
AREA=area(s)*1e6  # km to m2

writeRaster(AREA,paste('F:/Dataset/GC/area/AREA_1dot.tif',sep=""),options=c('TFW=YES'), overwrite=TRUE)

2.读取geotiff图像、计算全球排放的budget（from kg N ha^-1 yr^-1 to Tg N yr^-1）

R=georasterref('RasterSize',[360 720],'RasterInterpretation','cells',...
   'Latlim',[-90 90],'Lonlim',[-180 180],'ColumnsStartFrom','north');

Area=imread('F:/Dataset/GC/area/AREA_0dot5.tif');

%%这个是全球不同国家栅格图，1-204
str_WORLD='F:\Dataset\shp\worldtif\world_0dot5.tif';   
WORLD=double(imread(str_WORLD));
WORLD(WORLD>=255)=nan;
mmCoun=max(max(WORLD))

%%  crop NH3 emi，from 0.083 deg to 0.5 deg
CropNH3=imread(['F:\program\R\Ndep_Yield\data\raster\' 'CropNH3emi' num2str(2010) '.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/' 'CropNH3emi' num2str(2010) '.tif'],CropNH3,R); 

%%  calc buget
CropNH3_bug=CropNH3.*Area*1e-7;
sum(sum(CropNH3_bug))
geotiffwrite(['F:/program/R/Farmsize/data/' 'CropNH3emi' num2str(2010) 'bug.tif'],CropNH3_bug,R); 

%% summary by different country，矩阵思维
WORLD2=WORLD;
for ii=0:mmCoun
if(sum(sum(WORLD==ii))>0)
WORLD2(WORLD==ii)=nansum(CropNH3_bug(WORLD==ii));
end
end
geotiffwrite(strcat('F:/program/R/Farmsize/data/','CropNH3emi',num2str(2010),'bug_country.tif'),(WORLD2),R);

此外，使用R语言也可实现区域、国家统计

library(plyr)
library(rgeos)
library(rgdal)
library(raster)
library(rasterVis)
library(csvread)
library(viridis)
library(scales)
library(spatialEco)
library(Rcpp)
library(fasterize)
library(sf)

#############   stat  ratio NHx/NOy
shp_world=readOGR('K:/Dataset/shp/wordmap.shp')
shp_world_df=as.data.frame(shp_world,xy=TRUE)
bug_NH3=raster(paste('data/ratioNHx2NOy',2010,'.tif',sep="")) 

#summary by country to csv
sum_bug_NH3=zonal.stats(shp_world, bug_NH3, stats="mean")##sum
colnames(sum_bug_NH3)="sum_bug_NH3"
shp_world_df2=cbind(shp_world_df,sum_bug_NH3)
write.csv(shp_world_df2,paste("data/raster/sum_ratioNHx2NOy",2010,".csv",sep=""))

#summary by country to tif
shp_world@data$sum_bug_NH3=sum_bug_NH3$sum_bug_NH3 ###should be double
ext <- extent(-180,180,-90,90)
r <- raster(ext, res=0.083333333) 

shp_world2=st_as_sf(shp_world) #change to sf
r <- fasterize(shp_world2, r, field="sum_bug_NH3")
writeRaster(r,paste('data/raster/',"ratioNHx2NOy",2010,'_country.tif',sep=""),
            options=c('TFW=YES'), overwrite=TRUE)

注意zonal.stats不够精确，与Arcgis统计结果差别大，R语言的优势是统计分析（csv数据分析、统计函数分析等等），对空间栅格数据分析不够精确

标签：shp,sum,library,geotiff,求取,Matlab,NH3,world,bug
来源： https://www.cnblogs.com/rockman/p/16367887.html