clear all
close all

latitude= 19.5:0.5:54;
longitude=-130.5:0.5:-60;

file='udelairtemp_le.dat';
fid=fopen(file); 
var=fread(fid,142*70*348,'float32','ieee-le'); 
aa=find(abs(var)==32767);
var(aa)=NaN;
airtemp=reshape(var,142,70,348);
fclose(fid);

[ny,nx,ntime]=size(airtemp);

% read ascii file: load filename -ascii
%% mean and std

tempmean=mean(airtemp,3);
tempstd=std(airtemp,0,3);

%imagesc(flipud(airtemp(:,:,1)'));
%pcolor(airtemp(:,:,1)');
%% looping

for jj=1:ny
    for ii=1:nx
        for kk=1:ntime
            anomtemp(jj,ii,kk)=airtemp(jj,ii,kk)-tempmean(jj,ii);

        end
    end
end
%% show extreme values only
anomtemp2=zeros(ny,nx,ntime);
for jj=1:ny
    for ii=1:nx
        for kk=1:ntime
         
            if ( anomtemp(jj,ii,kk) <0) & (anomtemp(jj,ii,kk) >-10); 
%                  anomtemp2(jj,ii,kk)=NaN;
%              else if (anomtemp(jj,ii,kk) < 0); 
                anomtemp2(jj,ii,kk)=NaN;
            else
                anomtemp2(jj,ii,kk)=anomtemp(jj,ii,kk);
                %end
            end
        end
    end
end
%% save file
%save to binary format
file='sample.dat';
fid=fopen(file,'w');
count=fwrite(fid,airtemp,'float32');
fclose(fid);

%save to matlab format
save('sample.mat','airtemp');

%save to ascii format, use printf
%%
%fread(fileID,file size(in a string), precision, machine format)
%precision:specify the data form and size reading from fid
%machine format: specifies the order for reading bytes within the file. 
%after reading the variable in a string, you need to reshape the variable
%into a matrix format, lon*lat*time
%standard diviation flag: flag=0, normalize by (n-1); flag=1, normalize by n. n is the number of elements in the sample. 

%% plotting
figure(1)
axesm('MapProjection','mercator','FLatLimit',[20 54],...
  'FlonLimit',[-130 -60],'MLineLocation',20,...
  'PLineLocation',20,'MeridianLabel','on','ParallelLabel','on',...
  'LabelFormat','signed')
grid on;
framem on;

load coast
h = patchm(lat,long,'b','FaceColor','none');
%surfm(latitude,longitude,airtemp(:,:,13)');
contourm(latitude,longitude,airtemp(:,:,13)','color','black');
surfm(latitude,longitude,anomtemp2(:,:,1)');
states=shaperead('usastatehi','UseGeoCoords',true);
geoshow(states,'DefaultFaceColor','none','DefaultEdgeColor','black')
colorbar
%caxis([-10 40]);
title('Monthly air temperature (UDEL)','FontSize',16);
%print(1,'-djpeg','test.jpeg');
%%
%patchm: Project patch objects on map axes. 
%project and display patch (polygon) objects defined by their vertices given in lat and lon on the current map axes. 
%surfm: Project geolocated data grid on map axes