script demo: March 2010

Wednesday, March 31, 2010

Show Calendar

Tuesday, March 23, 2010

matlab: shaded world map based on coast data

function FillCoast(Is360,FillColor)
% To show a shaded worldmap based on the coast data coming with matlab
% Not perfectly, but does work some way.
% Usage:
%       FillCoast(Is360,FillColor)
% Copyright: http://scriptdemo.blogspot.com
if (nargin<1)
     Is360=0;
     FillColor='r';
elseif   (nargin==1)
     if (isnumeric(Is360) && length(Is360)==1)
        FillColor='r';
     else
        FillColor=Is360;
        Is360=0;
     end
elseif   (nargin>2)
     help FillCoast
     return
end
load coast;
if   ~ishold  hold on; end
if   Is360
     long(long<0)=long(long<0)+360;
     lat1=lat;
     lat1(long<0.5)=nan;
     plot(long,lat1,'k'); clear lat1
end
IndNaN=find(isnan(long));
Ns=1;
for   N=1:length(IndNaN)-1
      Ne=IndNaN(N)-1;
      lonseg=long(Ns:Ne);
      latseg=lat(Ns:Ne);
      Ns=Ne+2;
     if (length(lonseg)>3)
        N360=find(lonseg>340);
        NZero=find(lonseg<10);
        switch ~isempty(N360) + ~isempty(NZero)
           case {0,1}
                if (~Is360 && isempty(find(latseg>-60)))
                   % Antarctic
                   lonseg=[lonseg; lonseg(end); lonseg(1)];
                   latseg=[latseg;         -90;       -90];
                end
                hf=patch(lonseg,latseg,FillColor); set(hf,'linestyle','none')
           case {2}
                if isempty(find(latseg>-60))
                % Antarctic
                  Nmin=find(lonseg==min(lonseg));
                  lonseg=[lonseg(1:Nmin-1); lonseg(Nmin-1); lonseg(Nmin); lonseg(Nmin:end)];
                  latseg=[latseg(1:Nmin-1);         -90;       -90; latseg(Nmin:end)];
                  hf=patch(lonseg,latseg,FillColor); set(hf,'linestyle','none')
               else
                 lonseg1=lonseg(lonseg>200 & lonseg<=360);
                 if ~isempty(lonseg1)
                    lonseg1=[lonseg1;lonseg1(1)];
                    latseg1=latseg(lonseg>200 & lonseg<=360);
                    latseg1=[latseg1;latseg1(1)];
                    hf=patch(lonseg1,latseg1,'b'); set(hf,'linestyle','none')
                 end;
                 lonseg2=lonseg(lonseg<195);
                 if ~isempty(lonseg2)
                    lonseg2=[lonseg2;lonseg2(1)];
                    latseg2=latseg(lonseg<195);  latseg2=[latseg2;latseg2(1)];
                    hf=patch(lonseg2,latseg2,'g'); set(hf,'linestyle','none')
                 end
               end
           otherwise
             disp('You should not see this info...STOP')
        end
     end
end

Friday, March 12, 2010

NCO example III

#!/bin/bash
#  add a new row for the north pole with an average value of the last row
#  Usage:
#        NorthPole.sh ori_file.nc [ncvar mod_file.nc]
#        NetCdf Operator (NCO) package is required
#
GuessVar(){
if [ -e $1 ]; then
      ncvar=`ncdump -h $1 | grep "(time, lat, lon)" | awk '{print $2}' | awk -F\( '{print $1}'`
      if [ ${#ncvar} -eq 0 ]; then
        echo "Failed in guessing the variable name"
          exit
      else
          echo "   I guess the potential variable in the $1 is $ncvar"
      fi
else
      echo "$1 can not be found!!!"
      exit
fi
}
CheckLength()
{
NL=`echo $1 | awk -F\. '{print NF}'`   return $NL
}
#------------------------------------------------------------------------------------------------
if [ $# -lt 1 ]; then
    sed -n '2,7p' NorthPole.sh
    exit
else
   if [ $# -eq 1 ]; then
        if [ -e $1 ]; then
           ncfile_ori=$1
          ncfile_mod="mod_"${ncfile_ori}
          GuessVar $ncfile_ori
        else
           echo "$1 can not be found!!!"
           exit
       fi
     elif [ $# -eq 2 ]; then
       ncfile_ori=$1
      CheckLength $2
        if [ $? -eq 1 ]; then
         ncvar=$2
       ncfile_mod="mod_"${ncfile_ori}
       else
        ncfile_mod=$2
        GuessVar $ncfile_ori
      fi
   elif [ $# -eq 3 ]; then
        ncfile_ori=$1
        ncvar=$2
        ncfile_mod=$3
    fi
fi
echo "original nc file: $ncfile_ori, output ncfile: $ncfile_mod, variable: $ncvar"
myaddstr=" -s 'defdim(\"new_lat\",\$lat.size+1)' \
         -s 'new_var[\$time,\$new_lat,\$lon]=0.0f'\
       -s 'new_lat[\$new_lat]=0.0d'\
      -s 'new_lat(0:\$lat.size-1)=lat'\
          -s 'new_lat(\$lat.size)=90.0d' \
      -s 'new_var(:,0:\$lat.size-1,0:\$lon.size-1)=${ncvar}' \
      -s '*nnt=\$time.size; *nny=\$lat.size-1; *nnx=\$lon.size-1; *indt=0;' \
         -s 'for (indt==0; indt<nnt; indt++)  new_var(indt,nny+1,:)=${ncvar}(indt,nny,0:nnx).avg();'"
eval "ncap2 -O -v  ${myaddstr} ${ncfile_ori} ${ncfile_mod}"

#rename variables
ncwa -O -a lat ${ncfile_mod} tmp_${ncfile_mod}
ncks -O -x -v lat tmp_${ncfile_mod} ${ncfile_mod}
ncrename -d new_lat,lat -v new_var,${ncvar} -v new_lat,lat ${ncfile_mod} && rm -f tmp_${ncfile_mod}
ncatted -a standard_name,lat,o,c,"latitude" ${ncfile_mod}
ncatted -a long_name,lat,o,c,"latitude" ${ncfile_mod}
ncatted -a units,lat,o,c,"degrees_north" ${ncfile_mod}
ncatted -a axis,lat,o,c,"Y" ${ncfile_mod}

Thursday, March 11, 2010

NCO example II ncap2

addpole.nco:
defdim("new_y",$y.size+1);
new_uas[$new_y,$x]=0.0f;
new_y[$new_y]=0.0f;
new_y(0:$y.size-1)=y;
new_y($y.size)=90.0f;
new_uas(0:$y.size-1,0:$x.size-1)=uas;
*idx=0;
*meanpole=0.0;
while (idx<$x.size) {
meanpole+=uas($new_y.size-2,idx);
idx++;
}
meanpole=meanpole/idx;
new_uas($new_y.size-1,0:$x.size-1)=meanpole;

ncap2 -v -S addpole.nco ori.nc modified.nc