How can I export a matlab contour as a polygon shapefile?

I have a rainfall map for which I have created contour lines using the 'contour' function. I need to export the output contours as a shapefile.

Answer 1

[Walter Roberson]

This turned out to be easy.

First, the driver function:

function c2s_driver
    YourData = imread('cameraman.tif');

      cmatrix = contour(YourData);
      shapedata = contour2shape(cmatrix);

      shapewrite(shapedata, 'camera_contours.shp');
  end

and the workhorse:

function shp = contour2shape(cmatrix)
    %Converts a contour matrix to a shape structure
    %needs https://www.mathworks.com/matlabcentral/fileexchange/43162-c2xyz-contour-matrix-to-coordinates
    [x, y, z] = C2xyz(cmatrix);
    shp = struct('Geometry', 'PolyLine', 'X', x, 'Y', y, 'Z', num2cell(z));
end

Notice this needs C2xyz which you can download from the File Exchange or install using Add-On Explorer.

[Cedric]

This is brilliant! The only problems that I anticipate is that it won't manage the orientation correctly, which could be a problem when converting from polyline to polygon (if necessary).

[Walter Roberson]

I have no idea what the distinction is between Polygon and Polyline and PolygonM for this purpose ;-( Guess I should read https://www.esri.com/library/whitepapers/pdfs/shapefile.pdf

[Cedric]

Until we know the purpose, there is no need to solve this issue. If it is just for displaying lines in the GIS, your solution is worth 10 votes because it is simple.

[Anna Weeks]

Agreed that this is brilliant! Any ideas on the polyline to polygon conversion would be even more brilliant.

[Walter Roberson]

I think the below should create polygons properly:

function shp = contour2shape(cmatrix)
    %Converts a contour matrix to a shape structure
    %needs https://www.mathworks.com/matlabcentral/fileexchange/43162-c2xyz-contour-matrix-to-coordinates

      FixContourOrder = @(C) cellfun(@(v) [v(1), v(end:-1:1)], C, 'uniform', 0);

      [x, y, z] = C2xyz(cmatrix);
      shp = struct('Geometry', 'PolyGon', 'X', FixContourOrder(x), 'Y', FixContourOrder(y), 'Z', num2cell(z));
  end

The contour point order appears to be counter-clockwise, and the last point is not always the same as the first, so this fix-up reversed the point order and adds in the needed duplicate point.

I do not have a GIS to test the output with.

[Cedric]

I don't have much time sadly but .. little typo (corrected): in the call to STRUCT, it is FixContourOrder, and this cannot work directly unless Anna needs lines or polygons for display only and not for geoprocessing.

Let me first copy what I had started writing yesterday before I saw your post:

----------------------------------------------------------------------

Unless there is some function for this in the Mapping toolbox or in the FEX, this is not going to be completely trivial. If I had little time to find a solution and no access to the Mapping toolbox doc (which is the case because I am about to go to bed ;) and I am giving myself 10 minutes for writing a few hints), I would:

1. SHAPEREAD a random shapefile to get a geo struct, as a reminder of what I have to build for SHAPEWRITE

 >> S = shaperead( 'SomePolygons.shp' )
 S = 
  433×1 struct array with fields:
    Geometry
    BoundingBox
    X
    Y

 >> S(1)
 ans = 
  struct with fields:
       Geometry: 'Polygon'
    BoundingBox: [2×2 double]
              X: [-7.7717e+06 -7.7708e+06 -7.7756e+06 -7.7717e+06 NaN]
              Y: [3.0916e+06 3.0916e+06 3.0865e+06 3.0916e+06 NaN]

and remember that multi-part polygons are "NaN-separated". Also, the field BoundingBox is a 2x2 array that contains min/max X/Y:

 >> S(1).BoundingBox
 ans =
   1.0e+06 *
   -7.7756    3.0865
   -7.7708    3.0916

So if I get the polygons nodes coordinates, I can easily build BoundingBox. Seeing this, it seems doable to build polygons based on contours.

2. See whether I am able to access the contour data. I am building a test contour on a WGS1984 extent with 5 levels to see:

 >> [X,Y,Z] = peaks;
 >> figure
 >> [cContour,hContour] = contour( X*60, Y*30, Z, 5 ) ;

and looking at the contour object:

 >> hContour
 hContour = 
  Contour with properties:
    LineColor: 'flat'
    LineStyle: '-'
    LineWidth: 0.5000
         Fill: 'off'
    LevelList: [-4.1097 -1.6727 0.7643 3.2012 5.6382]
        XData: [49×49 double]
        YData: [49×49 double]
        ZData: [49×49 double]

Show all properties

            Annotation: [1×1 matlab.graphics.eventdata.Annotation]
          BeingDeleted: 'off'
            BusyAction: 'queue'
         ButtonDownFcn: ''
              Children: [0×0 GraphicsPlaceholder]
              Clipping: 'on'
         ContourMatrix: [2×432 double]
             CreateFcn: ''
                   ...

it seems that all we have to work with is in fact the ContourMatrix property, which happens to be the first output of CONTOUR, stored above in cContour (EDIT: maybe using CONTOURM, we would have access to special properties/method or functions of the Mapping tbx for exporting to shapefile). Plotting it we get:

which is not very encouraging, but we can use FEX by Kesh Ikuma or FEX by our colleague Chad Green to obtain nice X/Y coordinates for the contours.

3. Create polygons out of these contours, and here we have a problem, because most of these polygons are rings (like the multiple rings buffers from ArcGIS in a sense), which have an interior boundary. Now we can build these rings based on contiguous contours and detect the nesting using e.g. INPOLYGON, and then read the doc to know how to code interior rings (whether it is clockwise or counter-clockwise), but there is a step in complexity here.

Going to bed, but if I was you I would wait and see if someone comes with a solution already implemented. Otherwise it may be ~half a day of work.

[Cedric]

Ok, now here is why this cannot work directly. If you create polygons using contours, you don't create holes (interior boundaries), which means that polygons are not disjoint. We can see this importing a contour simpler than the cameraman:

 [X,Y,Z] = peaks;
 figure
 cmatrix = contour( X*60, Y*30, Z, 5 ) ;
 shapedata = contour2shape(cmatrix) ;
 shapewrite(shapedata, 'peaks_polygons.shp') ;

Now if we import that in a GIS, and use Z for defining polygons color, we get:

which seems fine, but if we set no color for the top polygon, we see that it is not:

The polygon below the dark one has no hole. This can not work for a series of geoprocessing operations (zonal stat, intersects, etc). There are ways to clean it though, nothing automatic as far as I remember unless you want to involve a little bit of Python (but look it up, don't take my work for that because I am programming most of my GIS processing and I am not a good user of the ArcGIS interface .. if it is what you are using).

Then if you keep exporting polylines ('PolyLine' in the call to SHAPEWRITE), you can convert polylines to polygons (in ArcGIS this is under "Data Management Tools/Features/Feature to Polygon") and you get the following, where I removed the filling of the top polygon:

This works, there is no overlay, but we loose Z or the attributes in the process.

[Walter Roberson]

Thanks, I fixed the typo.

I did not try to code insides with clockwise / counter-clockwise pairs. It might be a bit involved to figure out which polygon borders which.

I wonder if it would make more sense to use the z as a "measure" PolygonM type rather than an attribute? Or perhaps make it 3D points with constant Z ?

[Aleksey Tipikin]

What about polygons in the corner? They don't closes correctly after extraction. Here is the example.

%Latitude and longitude vectors
lat = -70:70;
lon = -180:180;
%full grid coordinates
[lat1,lon1] = ndgrid(lat,lon);
Z = sind(lat1*2)+sind(lon1*2);
%Georeference object
R=georefpostings([lat(1),lat(end)],[lon(1),lon(end)],size(lat1));
%Visualize maps
figure;
Map1 = worldmap('World');
setm(Map1,'FontSize',8,'MapProjection','mercator','MapLatLimit',[-70 70])
tightmap;
%Show Z
%contourfm(lat1,lon1,Z,25,'LineStyle','none');
contourfm(Z,R,25,'LineStyle','none');
contourcbar;
%Show coastline
geoshow('landareas.shp','FaceColor','none');
%Show specified level
[C,h] = contourfm(Z,R,[0.5 0.5]);
h.Children(3).FaceAlpha = 0;
h.Children(2).FaceAlpha = 0.4;
h.Children(2).FaceColor = 'w';
%Make polygon with coordinates of contour at specified level
[x, y, ~] = C2xyz(C);
g2 = geoshape(y{1},x{1},'Geometry','polygon');
geoshow(g2,'EdgeColor','k','FaceColor','r','FaceAlpha',0.4)