[Matlab] create a time series using datenum

function myT=yymmdd2x(inYear,inMonth,inDay) % create the time axis based on input time [year, mon, date], using datenum % usage: %           xTime=yymmdd2x(years, months, days); % or %           xTime=yymmdd2x(yyyymmdd); %  http://scriptdemo.blogspot.com if (nargin==0 || nargin>3)    help yymmdd2x;    return end if nargin==1    %yyyymmdd case    yyyymmdd=inYear;    inYear=floor(yyyymmdd/10000);    inDay=mod(yyyymmdd,100);    inMonth=floor(mod(yyyymmdd,10000)/100); clear yyyymmdd    totalDaysInYear=datenum(inYear+1,1,1)-datenum(inYear,1,1);    daysInYear=datenum(inYear,inMonth,inDay)-datenum(inYear,1,1)+1;    myT=daysInYear./totalDaysInYear+inYear; elseif nargin==2    % no date, set to 15th, may consider create 12 month for each year if necessary    if length(inYear)~=length(inMonth) && length(inMonth)<=12       % repeat inMonth for each year       inMonth=reshape(inMonth,1,[]);       mm=repmat(inMonth,1,length(inYear));       yy=reshape(repmat(reshape(inYear,1,[]),length(inMonth),1),1,[]);       clear inMonth inYear       totalDaysInYear=datenum(yy+1,1,1)-datenum(yy,1,1);       daysInYearA=datenum(yy,mm+1,1)-datenum(yy,1,1)+1;       daysInYear=datenum(yy,mm,1)-datenum(yy,1,1)+1;       myT=0.5*(daysInYearA+daysInYear)./totalDaysInYear+yy;    elseif (length(inYear)==length(inMonth))       totalDaysInYear=datenum(inYear+1,1,1)-datenum(inYear,1,1);       daysInYearA=datenum(inYear,inMonth+1,1)-datenum(inYear,1,1)+1;       daysInYear=datenum(inYear,inMonth,1)-datenum(inYear,1,1)+1;       myT=0.5*(daysInYearA+daysInYear)./totalDaysInYear+inYear;    else        disp('Too difficult for me to guess what your wanna to do...')        return    end else    totalDaysInYear=datenum(inYear+1,1,1)-datenum(inYear,1,1);    daysInYear=datenum(inYear,inMonth,inDay)-datenum(inYear,1,1)+1;    myT=daysInYear./totalDaysInYear+inYear; end

[Matlab] format a string to a specific length

function outStr=formatStr(inStr,nLen,opt) % pad the input string to specified length with blank space % usage: %           outStr=formatStr(inStr,expected-Length [,align-option]) %           if align-option == 0 'right-align' [ by default ] %                               == 1 'center-align' %                               == 2 'left-align' % history: %            Feb 24 2012: first coded for myShowNc %  http://scriptdemo.blogspot.com if nargin==1    nLen=0; elseif nargin==2    opt=0; elseif nargin~=3    help formatStr;    return end oriLen=length(inStr); if (oriLen>=nLen)     outStr=inStr; else     switch (opt)         case {0}             nSpace=nLen-oriLen;             outStr=[repmat(' ',1,nSpace),inStr];         case {1}           nSpace=ceil(0.5*(nLen-oriLen));           outStr=[repmat(' ',1,nSpace),inStr];           nSpace=nLen-oriLen-nSpace;           if (nSpace>0)              outStr=[outStr,repmat(' ',1,nSpace)];           end         case {2}             nSpace=nLen-oriLen;             outStr=[inStr,repmat(' ',1,nSpace)];         otherwise             disp('Unknown option');             return     end end

[Matlab] show the variable information from a netcdf file

function myShowNc(ncFileName) % show information of a given netcdf file, similar to ncdump %          using matlab native netcdf functions % usage: %           myShowNc(netcdf-filename) % http://scriptdemo.blogspot.com if nargin==0    help myShowNc    return end if ~exist(ncFileName,'file')    error([ncFileName,' is not found!']); end ncfid=netcdf.open(ncFileName,'NC_NOWRITE'); % check dimensions [ndims,nvars,ngatts,unlimdimID] = netcdf.inq(ncfid); myncinfo.unLimID=unlimdimID; myncinfo.dimname=cell(ndims,1); myncinfo.dimlength=zeros(ndims,1); myncinfo.varname=cell(nvars,1); myncinfo.vardims=cell(nvars,1); myncinfo.varscale=ones(nvars,1); myncinfo.varoffset=zeros(nvars,1); disp(['netcdf file: ',ncFileName]) disp(' dimension:') for Ndim=1:ndims     [myncinfo.dimname{Ndim},myncinfo.dimlength(Ndim)]=netcdf.inqDim(ncfid,Ndim-1);     disp([' ',fixStrLen(myncinfo.dimname{Ndim},10),': ',num2str(myncinfo.dimlength(Ndim))]) end disp(' variables:') disp([' ',fixStrLen('var-name |',12,0),fixStrLen('var-dims',30,0),fixStrLen('scale',10,1),fixStrLen('offset',10,1)]) % check variables for nv=1:nvars     [myncinfo.varname{nv},myncinfo.vardims{nv}, ...      myncinfo.varscale(nv),myncinfo.varoffset(nv)]=getVarInfo(ncfid,nv-1,myncinfo);     disp([' ', fixStrLen([myncinfo.varname{nv},' |'],12) , fixStrLen(myncinfo.vardims{nv},30), ...                 fixStrLen(num2str(myncinfo.varscale(nv)),10,1),fixStrLen(num2str(myncinfo.varoffset(nv)),10,1)]) end netcdf.close(ncfid); %-------------------------------------------------------------------------------------- function outStr=fixStrLen(inStr,nLen,opt) % renamed to formatStr in another post. % if  opt == 0  'right-align' [default case] %     opt == 1  'center-align' %     opt == 2  'left-align' if nargin==1    nLen=0; elseif nargin==2    opt=0; elseif nargin~=3    help fixStrLen;    return end oriLen=length(inStr); if (oriLen>=nLen)     outStr=inStr; else     switch (opt)         case {0}             nSpace=nLen-oriLen;             outStr=[repmat(' ',1,nSpace),inStr];         case {1}           nSpace=ceil(0.5*(nLen-oriLen));           outStr=[repmat(' ',1,nSpace),inStr];           nSpace=nLen-oriLen-nSpace;           if (nSpace>0)              outStr=[outStr,repmat(' ',1,nSpace)];           end         case {2}             nSpace=nLen-oriLen;             outStr=[inStr,repmat(' ',1,nSpace)];         otherwise             disp('Unknown option');             return     end end function [varname,dimsStr,myScale,myOffSet]=getVarInfo(ncfid,varID,myncinfo) [varname,vartype,vardimids,numAtts]=netcdf.inqVar(ncfid,varID); myScale=1; myOffSet=0; for numA=1:numAtts     tmpattname=netcdf.inqAttName(ncfid,varID,numA-1);     switch lower(tmpattname)         case {'scale_factor','scalefactor'}            myScale=netcdf.getAtt(ncfid,varID,tmpattname);         case {'add_offset','off_set','offset','offset_value'}            myOffSet=netcdf.getAtt(ncfid,varID,tmpattname);         otherwise            % do nothing      end end % get dimStr if ~isempty(vardimids)    dimsStr=myncinfo.dimname{vardimids(1)+1};    for nd=2:length(vardimids)          dimsStr=[myncinfo.dimname{vardimids(nd)+1},' x ',dimsStr];    end else    dimsStr='unlimited'; end

[Matlab] Convert a 2D mesh to a single kml file

function mesh2kml(kmlfile,long,lat,iskip,jskip) % convert a 2D mesh grid into a single plain kml file %             will create the horizontal and vertical mesh-line seperately % usage: %             mesh2kml(kml-filename,longitude,latitude [iskip, jskip]) %                     iskip: reduce the points when to generate i-direction-mesh-line %                     jskip: reduce the points when to generate j-direction-mesh-line %     e.g., %             lon=-90:2:20; %             lat=37:65; %             [lon,lat]=meshgrid(lon,lat); %             mesh2kml('mesh2kml_test',lon,lat,1,1); %             will produce a kml file named mesh2kml_test in current folder % %             http://scriptdemo.blogspot.com if nargin<3    help mesh2kml    return end if nargin==3    iskip=5;    jskip=5; elseif nargin==4    jskip=5; end fid=fopen([kmlfile,'.kml'],'wt'); writekmlHeader(fid,kmlfile); if length(long)~=numel(long) % mesh    [NY,NX]=size(long);    %i-index    longUD=flipud(long); latUD=flipud(lat);    longI=long; latI=lat;    longI(:,2:2:end)=longUD(:,2:2:end);    latI(:,2:2:end)=latUD(:,2:2:end); clear longUD latUD    jindex=sort(unique([1 jskip:jskip:NY NY]));    longI=reshape(longI(jindex,:),1,[]); latI=reshape(latI(jindex,:),1,[]);    fprintf(fid,'%s \n','<Placemark>');    fprintf(fid,'%s \n',['<name>',kmlfile,'-I</name>']);    fprintf(fid,'%s \n','<LineString>');    fprintf(fid,'%s \n','<tessellate>1</tessellate>'); %    fprintf(fid,'%s \n','<coordinates>');    fprintf(fid,'%0.6f,%0.6f,0.0 \n', [longI;latI]);    fprintf(fid,'%s \n','</coordinates>');    fprintf(fid,'%s \n','</LineString>');    fprintf(fid,'%s \n','</Placemark>');    %j-index    longLR=fliplr(long); latLR=fliplr(lat);    longJ=long; latJ=lat;    longJ(2:2:end,:)=longLR(2:2:end,:);    latJ(2:2:end,:)=latLR(2:2:end,:); clear longLR latLR    longJ=longJ'; latJ=latJ';    iindex=sort(unique([1 iskip:iskip:NX NX]));    longJ=reshape(longJ(iindex,:),1,[]); latJ=reshape(latJ(iindex,:),1,[]);    fprintf(fid,'%s \n','<Placemark>');    fprintf(fid,'%s \n',['<name>',kmlfile,'-J</name>']);    fprintf(fid,'%s \n','<LineString>');    fprintf(fid,'%s \n','<tessellate>1</tessellate>'); %    fprintf(fid,'%s \n','<coordinates>');    fprintf(fid,'%0.6f,%0.6f,0.0 \n', [longJ;latJ]);    fprintf(fid,'%s \n','</coordinates>');    fprintf(fid,'%s \n','</LineString>');    fprintf(fid,'%s \n','</Placemark>'); else    long=reshape(long,1,[]); lat=reshape(lat,1,[]);    fprintf(fid,'%s \n','<Placemark>');    fprintf(fid,'%s \n',['<name>',kmlfile,'</name>']);    fprintf(fid,'%s \n','<LineString>');    fprintf(fid,'%s \n','<tessellate>1</tessellate>'); %    fprintf(fid,'%s \n','<coordinates>');    fprintf(fid,'%0.6f,%0.6f,0.0 \n', [long;lat]);    fprintf(fid,'%s \n','</coordinates>');    fprintf(fid,'%s \n','</LineString>');    fprintf(fid,'%s \n','</Placemark>'); end writekmlEnd(fid); fclose(fid); function writekmlHeader(fid,kmlName) fprintf(fid,'%s \n','<?xml version="1.0" encoding="UTF-8"?>'); fprintf(fid,'%s','<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" '); fprintf(fid,'%s \n','xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">'); fprintf(fid,'%s \n','<Document>'); fprintf(fid,'%s \n',['<name>',kmlName,'.kml</name>']); function writekmlEnd(fid) fprintf(fid,'%s \n','</Document>'); fprintf(fid,'%s \n','</kml>');

[Matlab] Convert a line to a polygon or line kml file

function line2kml(kmlfile,long,lat,varargin) % convert a line to a kml line or polygon feature file with limited %              attributes support % usage: %              line2kml(kml-filename,longitude,latitude [, other opts]) %    opts: %              'style' : polygon or line %              'linecolor' : color of edge line %              'linealpha' : alpha value for edgecolor [1-255] %              'linewidth' : width of the edge line %              'fillcolor' : facecolor of the polygon %              'fillalpha' : alpha value for facecolor [1-255] %    e.g., %              np=1000; %              t=0:2*pi/np:2*pi; %              r=sin(t).*sqrt(abs(cos(t)))./(sin(t)+7/5)-2*sin(t)+2; %              long=r.*cos(t); lat=r.*sin(t); %              long=long/max(abs(long))*31-45; %              lat=lat/max(abs(lat))*35+38; %              line2kml('myheart',long,lat,'fillalpha',127,'style','polygon','linewidth',5,'linecolor','y','linealpha',255,'fillcolor','r'); %              will produce a kml file named myheart in current folder % %             @ http://scriptdemo.blogspot.com % demo case if nargin==0    np=1000;    t=0:2*pi/np:2*pi;    r=sin(t).*sqrt(abs(cos(t)))./(sin(t)+7/5)-2*sin(t)+2;    long=r.*cos(t); lat=r.*sin(t);    long=long/max(abs(long))*31-45;    lat=lat/max(abs(lat))*35+38;    line2kml('myheart',long,lat,'style','polygon','fillcolor','r','fillalpha',127,'linecolor','y','linewidth',5,'linealpha',255);    disp(['Please open the myheart.kml file and have a look in GoogleEarth!'])    return end if nargin<3    help line2kml    return end % check input arguments mykml.type='Polygon'; mykml.linecolor='r'; mykml.linecolorstr='r'; mykml.linewidth=1; mykml.fillcolor='r'; mykml.fillcolorstr='7f0000ff'; mykml.linealpha=127; mykml.fillalpha=127; while(size(varargin,2)>0)    switch lower(varargin{1})      case {'type','style'}           mykml.type=varargin{2};           varargin(1:2)=[];      case {'linewidth','lw'}           mykml.linewidth=varargin{2};           varargin(1:2)=[];      case {'linecolor','lc'}           mykml.linecolor=varargin{2};           varargin(1:2)=[];      case {'linealpha','linetransparency'}           mykml.linealpha=varargin{2};           varargin(1:2)=[];      case {'fillalpha','filltransparency','facealpha'}           mykml.fillalpha=varargin{2};           varargin(1:2)=[];      case {'fillcolor','fcolor'}           mykml.fillcolor=varargin{2};           varargin(1:2)=[];     otherwise           disp(['Unknown property: ',varargin{1}])           return    end end mykml.linecolorstr=getColor(mykml.linecolor,mykml.linealpha); mykml.fillcolorstr=getColor(mykml.fillcolor,mykml.fillalpha); fid=fopen([kmlfile,'.kml'],'wt'); writekmlHeader(fid,kmlfile); long=reshape(long,1,[]); lat=reshape(lat,1,[]); switch (mykml.type)   case {'Polygon','polygon','poly'}      fprintf(fid,'%s \n','<Style id="idpolyfill">');      fprintf(fid,'%s \n','<LineStyle>');      fprintf(fid,'%s \n',['<color>',mykml.linecolorstr,'</color>']);      fprintf(fid,'%s \n',['<width>',num2str(mykml.linewidth),'</width>']);      fprintf(fid,'%s \n','</LineStyle>');      fprintf(fid,'%s \n','<PolyStyle>');      fprintf(fid,'%s \n',['<color>',mykml.fillcolorstr,'</color>']);      fprintf(fid,'%s \n','</PolyStyle>');      fprintf(fid,'%s \n','</Style>');      fprintf(fid,'%s \n','<Placemark>');      fprintf(fid,'%s \n',['<name>',kmlfile,'</name>']);      fprintf(fid,'%s \n','<styleUrl>#idpolyfill</styleUrl>');      fprintf(fid,'%s \n','<Polygon>');      fprintf(fid,'%s \n','<tessellate>1</tessellate>'); %      fprintf(fid,'%s \n','<outerBoundaryIs>');      fprintf(fid,'%s \n','<LinearRing>');   case {'line'}      fprintf(fid,'%s \n','<Style id="idlineOnly">');      fprintf(fid,'%s \n','<LineStyle>');      fprintf(fid,'%s \n',['<color>',mykml.linecolorstr,'</color>']);      fprintf(fid,'%s \n',['<width>',num2str(mykml.linewidth),'</width>']);      fprintf(fid,'%s \n','</LineStyle>');      fprintf(fid,'%s \n','</Style>');      fprintf(fid,'%s \n','<Placemark>');      fprintf(fid,'%s \n',['<name>',kmlfile,'</name>']);      fprintf(fid,'%s \n','<styleUrl>#idlineOnly</styleUrl>');      fprintf(fid,'%s \n','<LineString>');      fprintf(fid,'%s \n','<tessellate>1</tessellate>'); %   otherwise    disp('Not defined') end fprintf(fid,'%s \n','<coordinates>'); %fprintf(fid,'%0.6f, %0.6f, 0.0 \n', [long;lat]); fprintf(fid,'%0.6f,%0.6f,0 \n', [long;lat]); % update for new version google earth, 2012.10 fprintf(fid,'%s \n','</coordinates>'); switch (mykml.type)   case {'Polygon','polygon','poly'}      fprintf(fid,'%s \n','</LinearRing>');      fprintf(fid,'%s \n','</outerBoundaryIs>');      fprintf(fid,'%s \n','</Polygon>');   case {'line'}      fprintf(fid,'%s \n','</LineString>');   otherwise    disp('Not defined') end fprintf(fid,'%s \n','</Placemark>'); writekmlEnd(fid); fclose(fid); function writekmlHeader(fid,kmlName) fprintf(fid,'%s \n','<?xml version="1.0" encoding="UTF-8"?>'); fprintf(fid,'%s','<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2" '); fprintf(fid,'%s \n','xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">'); fprintf(fid,'%s \n','<Document>'); fprintf(fid,'%s \n',['<name>',kmlName,'.kml</name>']); function writekmlEnd(fid) fprintf(fid,'%s \n','</Document>'); fprintf(fid,'%s \n','</kml>'); function colorHex=getColor(colorRGB,alpha) colorHex='00000000'; if ischar(colorRGB)    switch lower(colorRGB)      case {'r','red'}        colorHex(1:6)='FF0000';      case {'g','green'}        colorHex(1:6)='00FF00';      case {'b','blue'}        colorHex(1:6)='0000FF';      case {'k','black'}        colorHex(1:6)='000000';      case {'y','yellow'}        colorHex(1:6)='FFFF00';      case {'m','p'}        colorHex(1:6)='FF00FF';      otherwise           error(['Unkow color option: ',colorRGB])    end else    if (max(colorRGB)>1)       colorHex(1:2)=dec2hex(round(colorRGB(1)),2);       colorHex(3:4)=dec2hex(round(colorRGB(2)),2);       colorHex(5:6)=dec2hex(round(colorRGB(3)),2);    else       colorHex(1:2)=dec2hex(round(colorRGB(1)*255),2);       colorHex(3:4)=dec2hex(round(colorRGB(2)*255),2);       colorHex(5:6)=dec2hex(round(colorRGB(3)*255),2);    end end if alpha>1    colorHex(7:8)=dec2hex(round(alpha),2); else    colorHex(7:8)=dec2hex(round(alpha*255),2); end colorHex=colorHex([7 8 5 6 3 4 1 2]);

The demo case [run the script without any argument]: