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

2022-06-12 14:04:18 阅读：168 来源： 互联网

标签：shp sum library geotiff 求取 Matlab NH3 world bug

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

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

ICode9

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