How do I compute the area enclosed by contourm?

Hi all :)

I would like to know how to compute the area enclosed by a contour line plotted using contourm.

I have a 361x361 matrix respectivly for LAT, LON, sicstart. (sicEASE is also 361x361)

sicstart is basicly matrix containing 0, 1 or NaN.

I wanted to use polyarea(x,y) but I am not sure how to use it since I have a matrix of 0 ,1 or NaN and no (x,y) coordinates...

I have attached a figure of the contours I want to compute the area. Here is the code for the blue contour.

%Find ice at start week
sicstart           = sicEASE;
index1             = find(sicEASE >= 0.15);  %index for ice
index2             = find(sicEASE < 0.15);   %index for no ice
sicstart(index1)   = 1;                      %1 for ice
sicstart(index2)   = 0;                      %0 for no ice

figure(1)
ncpolarm('lat',66); hold on;
contourm(LAT,LON,sicstart,'b','LineWidth',1.5);  

Thanks for your help

Answer 1

[Star Strider]

I do not have the Mapping Toolbox, so I have no experience with contourm. However looking at the documentation, it appears to be the same as for contour, with the ‘C’ matrix in contourm being the same as the ‘M’ matrix in contour. (See the documentation on M for a guide to ite interpretation and to understand how it works.)

Recovering the disjointed contours is not trivial, however it can be done. An example is in how to fill the interior of a closed surface? . There are likely other examples available. This is the first one I found in my archive.

Then, use trapz or polyarea on the (x,y) vectors to calculate the areas of each, and sum them.

.

[A LL]

Hi Star Strider,

Thanks for your quick answer!

I am not familiar with the contour matrix but if I understood correclty I get the lat and lon of each enclosed contour in contourm.

Here is what I wrote and the new lines are commented as %new:

%Find ice at start week
sicstart           = sicEASE;
index1             = find(sicEASE >= 0.15);  %index for ice
index2             = find(sicEASE < 0.15);   %index for no ice
sicstart(index1)   = 1;                      %1 for ice
sicstart(index2)   = 0;                      %0 for no ice

figure(1)
ncpolarm('lat',66); hold on;
Cstart = contourm(LAT,LON,sicstart,'b','LineWidth',0.5);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Cstartlon       =   Cstart(1,:);                                %new
Cstartlat       =   Cstart(2,:);                                %new
Cstartlonfirst  =   Cstartlon(2:Cstart(2,1)+1);                 %new
Cstartlatfirst  =   Cstartlat(2:Cstart(2,1)+1);                 %new
afirst          =   polyarea(Cstartlonfirst,Cstartlatfirst);    %new

geoshow(Cstartlatfirst,Cstartlonfirst, 'DisplayType', 'point', 'Color', 'r');   %new

Here, Cstart(1,1) = 0.1000 and Cstart(2,1) = 3045

So, my first enclosed area would be defined by the coords given by the following entries of Cstart i.e. (Cstartlatfirst,Cstartlonfirst)

I attached the new figure I get.

I am still confused about this:

--> Can I still use polyarea for (lat,lon) coords instead of (x,y)?

--> Is there an efficient and simple way to use the contour matrix to get all the coords of the different contours? Since the contour matrix is divided in subset it seems tricky to do a simple loop...

Thanks

[Star Strider]

‘Here, Cstart(1,1) = 0.1000 and Cstart(2,1) = 3045’

That means that for the contour at 0.1000 there are 3045 (x,y) pairs for that particular segment.

‘Can I still use polyarea for (lat,lon) coords instead of (x,y)?’

Since I have no experience with the Mapping Toolbox functions, I cannot determine that. It will be necessary for you to experiment with it. Not enough information is provided for me to atempt to reconstruct the posted image.

The information in Create and Display Polygons appears to offer some options. From the links provided in that page, they appear to be compatible with the core MATLAB polygon functions, so polyarea could be an option.

‘Is there an efficient and simple way to use the contour matrix to get all the coords of the different contours? Since the contour matrix is divided in subset it seems tricky to do a simple loop.’

The loop is the only option I am aware of that will provide the (x.y) vectors for all the disjointed segments.

.

[A LL]

Hi Star Strider,

Thanks again for your help!

I ended up using a much simpler technique to compute my area.

The contour I plot is represented by a matrix of 0, 1 and NaN and I am only plotting the contour for the 1 values. I simply summed the number of grid cells containing the value 1 to get all the area covered by the contour. I also know the area of each grid cell so, it becomes super straightforward to get the area of the contour.

%find where there is ice on the EASE grid
sicstart           = sicEASE;
index1             = find(sicEASE >= 0.15);    %index for ice
index2             = find(sicEASE < 0.15);     %index for no ice
sicstart(index1)   = 1;                        %1 for ice
sicstart(index2)   = 0;                        %0 for no ice

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%new

%compute the total area of ice
area1cell   = 25*25;                     %[km^2] area of 1 grid cell on the EASE grid
area        = length(index1)*area1cell;  %sea ice area

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%plot figure
figure(1)
ncpolarm('lat',66); hold on;
contourm(LAT,LON,sicstart,'b','LineWidth',1.5);

Thanks again for your precious advice!

[Star Strider]

As always, my pleasure!

I am happy that you discovered as solution that works!