Asked by Anna Weeks
on 26 Sep 2017

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 by Walter Roberson
on 26 Sep 2017

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 Wannaz
on 26 Sep 2017

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 Wannaz
on 26 Sep 2017

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
on 26 Sep 2017

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 ?

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.