1 | function plot_gridded(md,data,options,plotlines,plotcols,i)
|
---|
2 | %PLOT_OVERLAY - superimpose radar image to a given field
|
---|
3 | %
|
---|
4 | % Usage:
|
---|
5 | % plot_gridded(md,options,plotlines,plotcols,i)
|
---|
6 | %
|
---|
7 | % See also: PLOTMODEL
|
---|
8 |
|
---|
9 | %process mesh and data
|
---|
10 | [x y z elements is2d isplanet]=processmesh(md,[],options);
|
---|
11 | [data datatype]=processdata(md,data,options);
|
---|
12 |
|
---|
13 | islevelset = exist(options,'levelset');
|
---|
14 | if islevelset
|
---|
15 | levelset = getfieldvalue(options,'levelset');
|
---|
16 | options2 = copy(options);
|
---|
17 | options2.removefield('caxis',false);
|
---|
18 | options2.removefield('log',false);
|
---|
19 | [levelset datatype]=processdata(md,levelset,options2);
|
---|
20 | end
|
---|
21 |
|
---|
22 | %check is2d
|
---|
23 | if ~is2d,
|
---|
24 | error('buildgridded error message: gridded not supported for 3d meshes, project on a layer');
|
---|
25 | end
|
---|
26 |
|
---|
27 | %Get xlim and ylim (used to extract radar image)
|
---|
28 | xlim=getfieldvalue(options,'xlim',[min(x) max(x)]);
|
---|
29 | ylim=getfieldvalue(options,'ylim',[min(y) max(y)]);
|
---|
30 |
|
---|
31 | isAxis = exist(options, 'axis');
|
---|
32 | if isAxis
|
---|
33 | myaxis = getfieldvalue(options,'axis');
|
---|
34 | if ~ischar(myaxis)
|
---|
35 | xlim = [myaxis(1), myaxis(2)];
|
---|
36 | ylim = [myaxis(3), myaxis(4)];
|
---|
37 | end
|
---|
38 | end
|
---|
39 |
|
---|
40 | postx=getfieldvalue(options,'posting',diff(xlim)/1000);
|
---|
41 | posty=getfieldvalue(options,'posting',diff(ylim)/1000);
|
---|
42 |
|
---|
43 | %Interpolating data on grid
|
---|
44 | x_m = xlim(1):postx:xlim(2);
|
---|
45 | y_m = ylim(1):posty:ylim(2);
|
---|
46 | data_grid=InterpFromMeshToGrid(elements,x,y,data,x_m,y_m,NaN);
|
---|
47 | data_grid_save = data_grid;
|
---|
48 | if size(data_grid,1)<3 | size(data_grid,2)<3,
|
---|
49 | error('data_grid size too small in plot_gridded, check posting and units');
|
---|
50 | end
|
---|
51 |
|
---|
52 | %Mask values if levelset>0
|
---|
53 | if islevelset
|
---|
54 | ls_grid=InterpFromMeshToGrid(elements,x,y,levelset,x_m,y_m,NaN);
|
---|
55 | data_grid(ls_grid>0) = NaN;
|
---|
56 | end
|
---|
57 |
|
---|
58 | %Process data_grid: add white in NaN and correct caxis accordingly
|
---|
59 | [data_nani data_nanj]=find(isnan(data_grid) | data_grid==-9999);
|
---|
60 | if exist(options,'caxis'),
|
---|
61 | caxis_opt=getfieldvalue(options,'caxis');
|
---|
62 | data_grid(find(data_grid<caxis_opt(1)))=caxis_opt(1);
|
---|
63 | data_grid(find(data_grid>caxis_opt(2)))=caxis_opt(2);
|
---|
64 | data_min=caxis_opt(1);
|
---|
65 | data_max=caxis_opt(2);
|
---|
66 | else
|
---|
67 | data_min=min(data_grid(:));
|
---|
68 | data_max=max(data_grid(:));
|
---|
69 | end
|
---|
70 |
|
---|
71 | %Select plot area
|
---|
72 | subplotmodel(plotlines,plotcols,i,options);
|
---|
73 |
|
---|
74 | %shading interp;
|
---|
75 | map = getcolormap(options);
|
---|
76 | image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(map)/(data_max-data_min))),map);
|
---|
77 | if exist(options,'shaded'),
|
---|
78 |
|
---|
79 | if exist(options,'dem'),
|
---|
80 | dem_grid=InterpFromMeshToGrid(elements,x,y,getfieldvalue(options,'dem'),x_m,y_m,NaN);
|
---|
81 | else
|
---|
82 | dem_grid=data_grid_save;
|
---|
83 | end
|
---|
84 | a = -45;
|
---|
85 | scut = 0.2;
|
---|
86 | c = 1;
|
---|
87 | % computes lighting from elevation gradient
|
---|
88 | [fx,fy] = gradient(dem_grid,x_m,y_m);
|
---|
89 | fxy = -fx*sind(a) - fy*cosd(a);
|
---|
90 | clear fx fy % free some memory...
|
---|
91 | fxy(isnan(fxy)) = 0;
|
---|
92 |
|
---|
93 | % computes maximum absolute gradient (median-style), normalizes, saturates and duplicates in 3-D matrix
|
---|
94 | r = repmat(max(min(fxy/nmedian(abs(fxy),1 - scut/100),1),-1),[1,1,3]);
|
---|
95 |
|
---|
96 | % applies contrast using exponent
|
---|
97 | rp = (1 - abs(r)).^c;
|
---|
98 | image_rgb = image_rgb.*rp;
|
---|
99 |
|
---|
100 | % lighter for positive gradient
|
---|
101 | k = find(r > 0);
|
---|
102 | image_rgb(k) = image_rgb(k) + (1 - rp(k));
|
---|
103 | end
|
---|
104 |
|
---|
105 | % set novalues / NaN to black color
|
---|
106 | if ~isempty(data_nani)
|
---|
107 | nancolor=getfieldvalue(options,'nancolor',[1 1 1]);
|
---|
108 | image_rgb(sub2ind(size(image_rgb),repmat(data_nani,1,3),repmat(data_nanj,1,3),repmat(1:3,size(data_nani,1),1))) = repmat(nancolor,size(data_nani,1),1);
|
---|
109 | end
|
---|
110 |
|
---|
111 | %plot grid
|
---|
112 | h=imagesc(xlim,ylim,image_rgb);
|
---|
113 | axis xy
|
---|
114 |
|
---|
115 | %last step: mesh gridded?
|
---|
116 | if exist(options,'edgecolor'),
|
---|
117 | A=elements(:,1); B=elements(:,2); C=elements(:,3);
|
---|
118 | patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
|
---|
119 | end
|
---|
120 |
|
---|
121 | %Apply options
|
---|
122 | if ~isnan(data_min) & ~isinf(data_min),
|
---|
123 | options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready
|
---|
124 | end
|
---|
125 | options=addfielddefault(options,'axis','xy equal'); % default axis
|
---|
126 | applyoptions(md,data,options);
|
---|
127 |
|
---|
128 | function y = nmedian(x,n)
|
---|
129 | %NMEDIAN Generalized median filter
|
---|
130 | % NMEDIAN(X,N) sorts elemets of X and returns N-th value (N normalized).
|
---|
131 | % So:
|
---|
132 | % N = 0 is minimum value
|
---|
133 | % N = 0.5 is median value
|
---|
134 | % N = 1 is maximum value
|
---|
135 |
|
---|
136 | if nargin < 2
|
---|
137 | n = 0.5;
|
---|
138 | end
|
---|
139 | y = sort(x(:));
|
---|
140 | y = interp1(sort(y),n*(length(y)-1) + 1);
|
---|