[7197] | 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
|
---|
[8472] | 10 | [x y z elements is2d isplanet]=processmesh(md,[],options);
|
---|
[8001] | 11 | [data datatype]=processdata(md,data,options);
|
---|
[7197] | 12 |
|
---|
[24762] | 13 | islevelset = exist(options,'levelset');
|
---|
| 14 | if islevelset
|
---|
| 15 | levelset = getfieldvalue(options,'levelset');
|
---|
[27402] | 16 | options2 = copy(options);
|
---|
[27356] | 17 | options2.removefield('caxis',false);
|
---|
| 18 | options2.removefield('log',false);
|
---|
[27402] | 19 | [levelset datatype]=processdata(md,levelset,options2);
|
---|
[24762] | 20 | end
|
---|
| 21 |
|
---|
[7197] | 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)]);
|
---|
[27540] | 30 |
|
---|
| 31 | isAxis = exist(options, 'axis');
|
---|
| 32 | if isAxis
|
---|
| 33 | myaxis = getfieldvalue(options,'axis');
|
---|
[27910] | 34 | if ~ischar(myaxis)
|
---|
| 35 | xlim = [myaxis(1), myaxis(2)];
|
---|
| 36 | ylim = [myaxis(3), myaxis(4)];
|
---|
| 37 | end
|
---|
[27540] | 38 | end
|
---|
| 39 |
|
---|
[20670] | 40 | postx=getfieldvalue(options,'posting',diff(xlim)/1000);
|
---|
| 41 | posty=getfieldvalue(options,'posting',diff(ylim)/1000);
|
---|
[7197] | 42 |
|
---|
| 43 | %Interpolating data on grid
|
---|
[21828] | 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);
|
---|
[15399] | 47 | data_grid_save = data_grid;
|
---|
[10413] | 48 | if size(data_grid,1)<3 | size(data_grid,2)<3,
|
---|
[14195] | 49 | error('data_grid size too small in plot_gridded, check posting and units');
|
---|
[10413] | 50 | end
|
---|
[7197] | 51 |
|
---|
[24762] | 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 |
|
---|
[7197] | 58 | %Process data_grid: add white in NaN and correct caxis accordingly
|
---|
[15399] | 59 | [data_nani data_nanj]=find(isnan(data_grid) | data_grid==-9999);
|
---|
[7197] | 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
|
---|
[12732] | 72 | subplotmodel(plotlines,plotcols,i,options);
|
---|
[7197] | 73 |
|
---|
| 74 | %shading interp;
|
---|
[15400] | 75 | map = getcolormap(options);
|
---|
[15399] | 76 | image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(map)/(data_max-data_min))),map);
|
---|
| 77 | if exist(options,'shaded'),
|
---|
[22701] | 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
|
---|
[15399] | 84 | a = -45;
|
---|
| 85 | scut = 0.2;
|
---|
| 86 | c = 1;
|
---|
| 87 | % computes lighting from elevation gradient
|
---|
[22701] | 88 | [fx,fy] = gradient(dem_grid,x_m,y_m);
|
---|
[15399] | 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));
|
---|
[14411] | 103 | end
|
---|
[15399] | 104 |
|
---|
| 105 | % set novalues / NaN to black color
|
---|
| 106 | if ~isempty(data_nani)
|
---|
[15400] | 107 | nancolor=getfieldvalue(options,'nancolor',[1 1 1]);
|
---|
[15399] | 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);
|
---|
[14411] | 113 | axis xy
|
---|
[7197] | 114 |
|
---|
| 115 | %last step: mesh gridded?
|
---|
| 116 | if exist(options,'edgecolor'),
|
---|
| 117 | A=elements(:,1); B=elements(:,2); C=elements(:,3);
|
---|
[11009] | 118 | patch('Faces',[A B C],'Vertices', [x y z],'FaceVertexCData',data_grid(1)*ones(size(x)),'FaceColor','none','EdgeColor',getfieldvalue(options,'edgecolor'));
|
---|
[7197] | 119 | end
|
---|
| 120 |
|
---|
| 121 | %Apply options
|
---|
[19026] | 122 | if ~isnan(data_min) & ~isinf(data_min),
|
---|
[15399] | 123 | options=changefieldvalue(options,'caxis',[data_min data_max]); % force caxis so that the colorbar is ready
|
---|
| 124 | end
|
---|
[15400] | 125 | options=addfielddefault(options,'axis','xy equal'); % default axis
|
---|
[7197] | 126 | applyoptions(md,data,options);
|
---|
[21240] | 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);
|
---|