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

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)

Answer 1

[Ameer Hamza]

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);

[Asanka Subasinghe]

Thank you for the reply.

I got an error message:

Error using min

Invalid option. Option must be 'omitnan' or 'includenan'.

Error in gwass1 (line 50)

contour_levels =[ min(Z,[],'all'):0.2:max(Z,[],'all')];

Aside from this, I was wondering if sorting the x and y coordinates misplaces the x-coordinates with their corresponding y values? Additionally, can you explain how the scattered interpolant function works? What exactly is it interpolating?

[Ameer Hamza]

It seems that you are using an older release of MATLAB (before R2019a). In that case, change this line to

contour_levels = min(Z,[],[1 2]):0.2:max(Z,[],[1 2]);

sorting x and y is just a requirement for making a meshgrid. It will not cause any misplacement of coordinates. scatteredInterpolant use the data scaledwellsx,scaledwellsy,gwtable and use it to give interpolated values on the meshgrid.

[Asanka Subasinghe]

Thank you for the explanation. Yes I am using R2017a.

I received another error.

Error using min

Dimension argument must be a positive integer scalar within indexing range.

Error in gwass1 (line 48)

contour_levels = min(Z,[],[1 2]):0.2:max(Z,[],[1 2]);

Does it have something to do with the []?

[Tommy]

Asanka,

Try this instead:

contour_levels = min(min(Z)):0.2:max(max(Z));

[Asanka Subasinghe]

Thanks Tommy, that seems to have worked.

[Asanka Subasinghe]

Is there a way for the contours to appear as if they are more spaced out? I want to display the contour labels but towards the top the contours are so closely spaced that the numbers do not appear very clearly as they are all bunched up.

Additionally, is there a way to show a flow direction of the wells?

Thank you

[Ameer Hamza]

Asanka, you can increase the distance between the contour lines by increasing the distance between their levels

contour_levels = min(min(Z)):0.5:max(max(Z)); % 0.5 distance between lines

[Tommy]

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));

[Asanka Subasinghe]

Thanks to the both of you for the help