How to create contour of water table at 0.2 m intervals?

5 次查看（过去 30 天）
Asanka Subasinghe
Asanka Subasinghe 2020-4-13
评论： Asanka Subasinghe 2020-4-14
I have some wells that are located by a lake. I have the ground watertable elevation at each well, as well as the elevation at the lake. I created x and y coordinates for these elevations and am now trying to combine it to create contour lines of the water table. How can I do this at 0.2 m intervals?
Not sure what the problem is with using the reshape function here as well.
clc;clear all;close all
% x and y coordinates of well locations and lake points in cm (3.4cm/km)
w1=[5.9,5.6];
w8=[8,7.6];
w17=[11.9,2];
w21=[9.7,5.3];
w27=[4.6,4.6];
w56=[5.9,1.7];
w58=[5.9,1];
w59=[6.7,4.1];
lakepoint1=[4,6.55];
lakepoint2=[5.75,6.9];
lakepoint3=[7.7,7.65];
% Combined coordinates into a single matrix
wells=[w1;w8;w17;w21;w27;w56;w58;w59;lakepoint1;lakepoint2;lakepoint3];
% Converted well locations to metres
scaledwells=wells*0.441176*1000;
%Groundwater table elevations at each well and at 3 points of the lake. (in metres)
gwtablew1=236.07;
gwtablew8=236.65;
gwtablew17=240.37;
gwtablew21=238.13;
gwtablew27=235.96;
gwtablew56=236.66;
gwtablew58=237.70;
gwtablew59=237.90;;
gwtablelake1=234.19;
gwtablelake2=234.19;
gwtablelake3=234.19;
%Combined into single matrix below
gwtable=[gwtablew1;gwtablew8;gwtablew17;gwtablew21;gwtablew27;gwtablew56;gwtablew58;gwtablew59;gwtablelake1;gwtablelake2;gwtablelake3];
%created scatter plot of x and y coordinates of wells and lake points)
scatter(scaledwells(:,1),scaledwells(:,2))
xlabel('m')
ylabel('m')
grid on
%Trying to produce contours
scaledwellsx=scaledwells(:,1)
scaledwellsy=scaledwells(:,2)
nx=length(scaledwellsx);
ny=length(scaledwellsy);
z=reshape(gwtable,ny,nx)
I am receiving the error:
Error using reshape
To RESHAPE the number of elements must not change.
Error in gwass1 (line 45)
z=reshape(gwtable,ny,nx)

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

采纳的回答

Ameer Hamza
Ameer Hamza 2020-4-13
You first need to interpolate your data to a grid before plotting contour. Try the following code
clc;clear all;close all
% x and y coordinates of well locations and lake points in cm (3.4cm/km)
w1=[5.9,5.6];
w8=[8,7.6];
w17=[11.9,2];
w21=[9.7,5.3];
w27=[4.6,4.6];
w56=[5.9,1.7];
w58=[5.9,1];
w59=[6.7,4.1];
lakepoint1=[4,6.55];
lakepoint2=[5.75,6.9];
lakepoint3=[7.7,7.65];
% Combined coordinates into a single matrix
wells=[w1;w8;w17;w21;w27;w56;w58;w59;lakepoint1;lakepoint2;lakepoint3];
% Converted well locations to metres
scaledwells=wells*0.441176*1000;
%Groundwater table elevations at each well and at 3 points of the lake. (in metres)
gwtablew1=236.07;
gwtablew8=236.65;
gwtablew17=240.37;
gwtablew21=238.13;
gwtablew27=235.96;
gwtablew56=236.66;
gwtablew58=237.70;
gwtablew59=237.90;;
gwtablelake1=234.19;
gwtablelake2=234.19;
gwtablelake3=234.19;
%Combined into single matrix below
gwtable=[gwtablew1;gwtablew8;gwtablew17;gwtablew21;gwtablew27;gwtablew56;gwtablew58;gwtablew59;gwtablelake1;gwtablelake2;gwtablelake3];
%created scatter plot of x and y coordinates of wells and lake points)
scatter(scaledwells(:,1),scaledwells(:,2))
xlabel('m')
ylabel('m')
grid on
hold on
%Trying to produce contours
scaledwellsx=scaledwells(:,1);
scaledwellsy=scaledwells(:,2);
[X,Y] = meshgrid(sort(scaledwellsx), sort(scaledwellsy));
F = scatteredInterpolant(scaledwellsx,scaledwellsy,gwtable);
Z = F(X,Y);
contour_levels = min(Z,[],'all'):0.2:max(Z,[],'all');
contour(X, Y, Z, contour_levels);
  9 个评论
显示 7更早的评论
Tommy
Tommy 2020-4-14
Alternatively, if you want to keep the contour lines at 0.2 m intervals as initially requested, you can use the clabel function to add your labels. clabel allows you to specify which contour lines should be labeled. For example, to label only every 4th contour line:
contour_levels = min(min(Z)):0.2:max(max(Z));
[M, c] = contour(X, Y, Z, contour_levels);
clabel(M, c, contour_levels(1:4:end));
or to label only the upper half of contour lines:
clabel(M,c, contour_levels(end/2:end));
You might have to play around with it to get the labels exactly where you want.
Below was my attempt to find and plot the flow directions using gradient:
[DX,DY] = gradient(Z);
zi = arrayfun(@(i) find(X(1,:)==scaledwellsx(i)&Y(:,1)==scaledwellsy(i),1), 1:numel(scaledwellsx));
% ^ linear indices within Z (and DX, DY) of each well/lake
quiver(scaledwellsx',scaledwellsy',-DX(zi),-DY(zi));

请先登录，再进行评论。

更多回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by