Grid to time and Grid to lon/lat conversion

6 次查看（过去 30 天）
Tanziha Mahjabin
Tanziha Mahjabin 2020-10-20
评论： dpb 2020-10-21
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')
  1 个评论
Walter Roberson
Walter Roberson 2020-10-20
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Tanziha Mahjabin
Tanziha Mahjabin 2020-10-21
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")
  3 个评论
显示 1更早的评论
dpb
dpb 2020-10-21
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Data Distribution Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by