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

一、Matlab区域统计

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_0dot5.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数据分析、统计函数分析等等），对空间栅格数据分析不够精确

extract与arcgis、matlab的方法计算结果基本一致，但是速度较慢，本方法也可使用（如对速度要求不高）

bug_NH3=raster(paste('F:/program/R/Ndep_sat/data/yr/Buget_NHx',yr,'.tif',sep="")) 
  #summary by country to csv
  sum_bug_NH3=extract(bug_NH3,shp_world, fun=sum)
  colnames(sum_bug_NH3)="sum_bug_NH3"
  shp_world_df2=cbind(shp_world_df,sum_bug_NH3)
  write.csv(shp_world_df2,paste("F:/program/R/Ndep_sat/data/yr/sum_bug_NHx",yr,".csv",sep=""))

三、用Arcpy区域统计，然后用Matlab或者R对dbf文件进项合并，分析等操作，强烈推荐！计算结果准确

import sys, string, os, arcgisscripting
gp = arcgisscripting.create()
gp.CheckOutExtension("spatial")
gp.AddToolbox(r"C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx")
gp.overwriteOutput=1
import arcpy
from arcpy.sa import *
from arcpy import env as myEnv

# Set environment settings
myEnv.workspace = r"F:/program/R/Ndep_sat/data"

inZoneData = r"F:/Dataset/shp/wordmap.shp"
zoneField = "NAME"

inValueRaster = r"F:/program/R/Ndep_sat/data/yr/Buget_NHx2010.tif"
outTable=str(r"F:/program/R/Ndep_sat/data/yr"+"/" + "sum_bug_NHx_2010.dbf")
outZSaT = ZonalStatisticsAsTable(inZoneData, zoneField, inValueRaster, outTable, "NODATA", "SUM")

outZSaT = ZonalStatistics(inZoneData, zoneField, inValueRaster, "SUM","NODATA",)
outZSaT.save(r"F:/program/R/Ndep_sat/data/yr"+"/" + "sum_bug_NHx_2010v2.tif")