Grid to time and Grid to lon/lat conversion

Hi,

I have rotated the current grid and then plot. But the plot is giving on grid for both axis. How can I show (lat,lon) on y axis and time on x axis.

I tried to use datetick; but it's giving wrong date.

(As the grid is rotated; everything is calculated acordingly. Also because of grid rotation there should be lat and lon both on y axis.)

[NTR NIR]=size(UTR);
[NTR NIR]=meshgrid(1:NIR,1:NTR);
figure(4); clf; quiver(NIR,NTR,UTR,VTR,1,'k')
% set(gca,'xlim',[ min(TIMEw) max(TIMEw)])
datetick('x','mm/dd/yy','keeplimits','keepticks')

[Walter Roberson]

It is not clear that it makes sense to display a time axes in that plot. Your X and Y are NIR and NTR which go up to about 320 and 55 respectively. Is the 320 intended to be days? If so then for mm/dd/yy to be meaningful there would need to be a base date to start from.

Answer 1

[Tanziha Mahjabin]

The full code was something like the following. I did rotate the grid and calculated al variables accordingly. Then the weekly average is done. I want to see the change of the current over time. So I need the time on x axis. The figure I uploaded before was showing 304 time steps.

% find a line orthogonal to the site-to-site baseline in its midpoint; then
% find the indexes of the radar grid points closest to this orthogonal line

% some useful tools
clear all; clc; close all
% addpath('~/Documents/Github/scosoli/tmp/acorn_wera_Radial_DM-mbp/Helpers/HFRadarmap4_1/'); 
path(path,'H:\Data\Radar_job\acorn_daily_averaging\')

fileIn='IMOS_aggregation_20200124T074252Z.nc';
time=ncread(fileIn, 'TIME');
reftime=ncreadatt(fileIn,'TIME','units'); reftime=reftime(12:end-4);
reftime=datenum(reftime,'yyyy-mm-dd HH:MM:SS');
time=time+reftime;
% [iw,jw] = ndgrid(iwk,1:size(sst,2));

                                                 
U=single(ncread(fileIn, 'UCUR'));
V=single(ncread(fileIn, 'VCUR'));

Usd=single(ncread(fileIn, 'UCUR_sd'));
Vsd=single(ncread(fileIn, 'VCUR_sd'));

Uqc=single(ncread(fileIn, 'UCUR_quality_control'));
% Vqc=ncread(fileIn, 'VCUR_quality_control');
                                                
% get variable size and prefilter based on available QC flags for 1-h data
[nJ nI nT]=size(U);
id=Uqc>2;
u=U; u(id)=nan;
v=V; v(id)=nan;

LON=ncread(fileIn,'LONGITUDE');
LAT=ncread(fileIn,'LATITUDE');

% get sizes
[nJ,nI,nT]=size(U);


% typically for WERAs we only store the unique lon-lat vectors rather than
% a grid
if min(size(LON))==1
    [LAT LON]=meshgrid(LAT,LON);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Y M D H MN S]=datevec(time);
% average values per day
% [ad,~,cd] = unique(data(:,1:3),'rows');
% out_day = [ad,accumarray(cd,data(:,5),[],@nanmean)];

[ad,~,cd] = unique([Y M D],'rows');
T=(datenum([ad(:,1) ad(:,2) ad(:,3) zeros(size(ad(:,3))) zeros(size(ad(:,3))) zeros(size(ad(:,3)))]));

%  preallocation of variables
TIMED=T; 
    UD=  nan(nJ,nI,size(ad,1));
    VD=  nan(nJ,nI,size(ad,1));
  UsdD=  nan(nJ,nI,size(ad,1)); 
  VsdD=  nan(nJ,nI,size(ad,1)); 
NOBS1D=  nan(nJ,nI,size(ad,1)); 
NOBS2D=  nan(nJ,nI,size(ad,1)); 
   QCD=  nan(nJ,nI,size(ad,1)); 

for j=1:nJ; 
    for i=1:nI;
        if ~all(isnan(squeeze(u(j,i,:))))
            clear out_U   ; out_U   = [ad,accumarray(cd,squeeze(u(j,i,:)),[],@nanmean )];
            clear out_Usd ; out_Usd = [ad,accumarray(cd,squeeze(u(j,i,:)),[],@nanstd  )];
            clear out_V   ; out_V   = [ad,accumarray(cd,squeeze(v(j,i,:)),[],@nanmean )];
            clear out_Vsd ; out_Vsd = [ad,accumarray(cd,squeeze(v(j,i,:)),[],@nanstd  )];
            clear out_N   ; out_N   = [ad,accumarray(cd,squeeze(u(j,i,:)),[],@nancount)];
    
              UD(j,i,:) = out_U(:,4)  ;
            UsdD(j,i,:) = out_Usd(:,4);
              VD(j,i,:) = out_V(:,4)  ;
            VsdD(j,i,:) = out_Vsd(:,4);
          NOBS1D(j,i,:) = out_N(:,4)  ;
          NOBS2D(j,i,:) = out_N(:,4)  ;

        end
    end
end



id=NOBS1D<=2            ; QCD(id)=4;
id=NOBS1D>2 & NOBS1D<= 6; QCD(id)=3;
id=NOBS1D>6 & NOBS1D<=12; QCD(id)=2;
id=NOBS1D>12            ; QCD(id)=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load('variables_bonc.mat')
% locations of the radar stations
radialSite(1).location=[139.8497833 -37.3285500]; %NOCR
radialSite(2).location=[140.4562167 -37.9392667]; %BFCV

%% second approach: remap the u,v currents on a different grid rotated and
% parallel to the baseline, then get a 'conventional' transect centered 
% at the lon-lat of the max data return

% new grid
[x12 y12] = lonlat2km(radialSite(1).location(1),radialSite(1).location(2),...
                      radialSite(2).location(1),radialSite(2).location(2));
d12 = sqrt(x12^2+y12^2);

 east = cos(atan(y12/x12))*d12/2;
north = sin(atan(y12/x12))*d12/2;
% questa coppia (lon,lat) definisce il centro della baseline tra sit1 e
% sit2. tutta la griglia viene generata relativamente a questo punto.

[lon, lat] = km2lonlat(radialSite(1).location(1),radialSite(1).location(2),...
                      east,north);

% now define the grid relatively to the centre of the baseline connecting
% site1 to site2. the grid will be parallel to the baseline.
grid1 = 1;
grid2 = 0;

if grid1
    grdRes = 6;
    xGrid = -34*6:grdRes:34*6; % questa nella configurazione originale e nelle
                           %  varianti -- ad eccezione dell'ultima
    yGrid = -34*6:grdRes:34*6;
end
if grid2
    grdRes = 1;
    xGrid = -16:grdRes:16; % questa nella configurazione originale e nelle
                           %  varianti -- ad eccezione dell'ultima
    yGrid = -20:grdRes:-1;
end

[x y] = meshgrid(xGrid,yGrid);

[xGr,yGr] = rot_alfa(x,y,-atand(y12/x12));% la rotazione della griglia 
                                          % riporta tutto parallelo alla
                                          % baseline sit1-sit2
% [lonGr,latGr]=km2lonlat(lon,lat,...
%                             xGr,yGr);

[lonGr,latGr]=km2lonlat(mean(LON(:)),mean(LAT(:)),...
                            xGr,yGr);

% load('variables_bonc.mat')

% regrid everything
UR=nan(size(UD));
VR=nan(size(UD));

for t=1:length(TIMED)
    UR(:,:,t)=griddata(LON,LAT,squeeze(UD(:,:,t)),lonGr,latGr);
    VR(:,:,t)=griddata(LON,LAT,squeeze(VD(:,:,t)),lonGr,latGr);
end
% find the point with the max data return
[nJR nIR nTR]=size(UR);
perGoodR=zeros(nJR,nIR);
for jR=1:nJR
    for iR=1:nIR
        if ~all(isnan(squeeze(UR(jR,iR,:))))
            perGoodR(jR,iR)=numel(find(~isnan(squeeze(UR(jR,iR,:)))))/nTR*100;
        end
    end
end

[v pIR]=max(max(perGoodR));
[v pJR]=max(max(perGoodR'));

figure(3); clf; hold on
pcolor(lonGr,latGr,perGoodR)
plot(lonGr(:,pIR),latGr(:,pIR),'gd');


clear uTR vTR
% for t=20:-1:1
for t=1:length(TIMED)
    uTR(t,:)=squeeze(UR(:,pIR,t));
    vTR(t,:)=squeeze(VR(:,pIR,t));
end


% now I have everything in daily basis data on rotated grid. I will convert
% the data on weekly basis now

% get variable size and prefilter based on available QC flags
[nJ nI nT]=size(uTR);
[nJ nI nT]=size(vTR);

% step 1: sort out time and week days
iwk = cumsum([1;diff(weekday(TIMED) == 2) == 1]);
% [iw,jw] = ndgrid(iwk,1:size(sst,2));
[iw,jw] = ndgrid(iwk,1:size(uTR,2));
[iw,jw] = ndgrid(iwk,1:size(vTR,2));

z = accumarray(iwk,TIMED,[],@(x){strjoin(cellstr(datestr(x([1,end]),...
                                                     'dd-mmm-yyyy')),' ... ')});

% step 2: calculate averages and get number of obs. for each average. need
% to use cell2mat as data is stored into MATLAB cell format
% NB num2cell is required as there is a string (char) with time information that
% needs to be managed. conversion to cell handles that
clear wkly; wkly = [z, num2cell(accumarray([iw(:),jw(:)],uTR(:),[],@nanmean))];
UTR=cell2mat(wkly(:,2:end));
clear wkly; wkly = [z, num2cell(accumarray([iw(:),jw(:)],vTR(:),[],@nanmean))];
VTR=cell2mat(wkly(:,2:end));


% get time
TIMEw=cell2mat(wkly(:,1));
% tvec=reshape(TIMEw,nweeks, nJ, nI)

TIMEw=datenum(TIMEw(:,1:11));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[NTR NIR]=size(UTR);
[NTR NIR]=meshgrid(1:NIR,1:NTR);

figure(4); clf; quiver(NIR,NTR,UTR,VTR,1,'k')

% set(gca,'xlim',[ min(TIMEw) max(TIMEw)])
datetick('x','mm/dd/yy','keeplimits','keepticks')

% lon_C2=118.44;
% lat_C2=-19.99;
% [lon_thin,lat_thin]=meshgrid(lon,lat);
% dC2=sqrt((lon_thin-lon_C2).^2+(lat_thin-lat_C2).^2);
% [ixC2,iyC2]=find(dC2==min(min(dC2)));
% 
% z=data.frame(
%   x=seq(as.Date("2010-07-16"), by="+1 month", length.out=10)
%   y=1:10
% )
% plot(y~x, data=z)
% abline(v=axis.Date(1, z$x), col="grey80")

[Walter Roberson]

We would need IMOS_aggregation_20200124T074252Z.nc to test with.

When I glance at the code (without being able to test it), it does not appear to me that you have representable time on your x axis.

[Tanziha Mahjabin]

Hi Walter,

I have attached the file in dropbox.

https://www.dropbox.com/s/us2h51zisl0pq6h/IMOS_aggregation_20200124T074252Z.nc?dl=0

[dpb]

time=ncread(fileIn, 'TIME');
reftime=ncreadatt(fileIn,'TIME','units'); reftime=reftime(12:end-4);
reftime=datenum(reftime,'yyyy-mm-dd HH:MM:SS');
time=time+reftime;

Don't use deprecated datenum, use datetime instead. Then plot routines are datetime aware and you'll get time axes automagically if using a datetime variable as the abscissa.

If the time is some number of time units from a reference, it will be a duration added to the reference datetime