clc
clear
ncFilePath='E:\批處理讀取NC檔案-MATLAB語言\1NC檔案\VIC_Output_China_evap_1952.nc';


%% 顯示結構
% ncdisp(ncFilePath);%顯示nc檔案的所有結構，以便大概瞭解裡面的內容
% ncdisp(ncFilePath,'evap');%顯示指定變數的內容，注意一定要是變數variables才可以
% ncdisp(ncFilePath,'/','min');%簡單顯示結構以及定義
% ncdisp(ncFilePath,'/','full');%全部顯示所有結構和定義資訊


%% 讀取變數值
% ncid = netcdf.open(ncFilePath,'NOWRITE'); %開啟nc檔案返回索引ID
% [ndims,nvars,ngglobalatts,unlimdimid] = netcdf.inq(ncid);%獲取維數，變數數，全域性屬性數量，
% [varname,xtype,dimids,natts] = netcdf.inqVar(ncid,0); %根據變數索引號獲取變數的名稱


lon=ncread(ncFilePath,'lon');%讀取經度變數
lat=ncread(ncFilePath,'lat');%讀取緯度變數
time=ncread(ncFilePath,'time');%讀取時間變數
evap=ncread(ncFilePath,'evap');%獲取蒸散發變數資料
evap1=evap(:,:,200);
evap1(evap1==0)=[];
aa=evap1;
[h,l]=find(aa==0);
aa(h,l)=[];
%% 顯示資料
% pcolor(lat,lon,evap1);
% [x,y]=meshgrid(lon,lat);%根據經緯度資訊產生格網，284列（經度），164列（緯度）
% phandle=pcolor(x,y,evap1');%顯示一個矩陣，其中x,y,evap1的行列數必須一致。類似surface函式
% colorbar
% imwrite(evap1','E:\批處理讀取NC檔案-MATLAB語言\aa.tif','tif')
%% 儲存為地理柵格tif格式，以便arcgis讀取
data=flipud(evap1');
R = georasterref('RasterSize', size(data),'Latlim', [double(min(lat)) double(max(lat))], 'Lonlim', [double(min(lon)) double(max(lon))]);%地理柵格資料參考物件(類) 
geotiffwrite('E:\批處理讀取NC檔案-MATLAB語言\aa4.tif',data,R);