All,
We try to generate contour polygons from an unstructured triangular grid
stored in a netcdf file:
import netCDF4
import matplotlib.tri as tri
# read data
var = netCDF4.Dataset('filename.cdf').variables
x = var['x'][:]
y = var['y'][:]
elems = var['element'][:,:]-1
data = var['attrname'][:]
# create triangulation
triang = tri.Triangulation(x, y, triangles=elems)
# generate contour polygons
levels = numpy.linspace(0, maxlevel, num=numlevels)
contour = plt.tricontourf(triang, data, levels=levels)
# extract geometries
for coll in contour.collections:
# handle all geometries for one level
for p in coll.get_paths():
polys = p.to_polygons()
...
At this point we assume, that polys[0] is a linear ring to be interpreted as
a polygon exterior and polys[1:] are the corresponding interiors for
polys[0].
Here are our questions:
Is this assumption correct?
Is there any detailed documentation describing the structure of the returned
geometries?
Are the linear rings supposed to be correctly oriented (area > 0 for
exteriors and area < 0 for the interiors)?
Thanks!
Regards Hartmut
---------------
http://boost-spirit.com
http://stellar.cct.lsu.edu
On 26 October 2014 00:18, Hartmut Kaiser <har...@gm...> wrote: > At this point we assume, that polys[0] is a linear ring to be interpreted > as > a polygon exterior and polys[1:] are the corresponding interiors for > polys[0]. > > Here are our questions: > > Is this assumption correct? > Is there any detailed documentation describing the structure of the > returned > geometries? > Are the linear rings supposed to be correctly oriented (area > 0 for > exteriors and area < 0 for the interiors)? > Hello Hartmut, In brief, the answers are no, no and yes. In more detail, assuming polys is not empty then it will contain one or more polygon exteriors and zero or more interiors, and they can be in any order. Here is a simple example where polys[0] is an interior and polys[1] an exterior: x = [0, 0, 1, 1, 0.5] y = [0, 1, 0, 1, 0.5] z = [0.5, 0.5, 0.5, 0.5, 0] triang = tri.Triangulation(x, y) contour = plt.tricontourf(triang, z, levels=[0.2, 0.4]) The returned geometries are purposefully not documented. They are an 'implementation detail' and not considered part of the public interface. and as such they could change at any time and hence should not be relied upon. Of course you can choose to access them if you wish, as I do myself sometimes, but we make no promises about what the order of the polygons is, or that it won't change tomorrow. In reality the order of the polygons is chosen to be something that is easy for both the contour generation routines to create and for the various backends to render. If you were to look at the output generated by contourf, you will see that it is organised differently from that produced by tricontourf and is more like you would like it to be, i.e. one or more groups of an exterior polygon followed by zero or more interiors. This is historical as the contourf code dates from before all of the backends were able to render arbitrary groups of exterior and interior polygons, and so the contourf code has to calculate the order for the backends. When the tricontourf code was written the backends were all able to calculate the exterior/interior containment themselves, so there was no need for tricontourf to do it as well. Ian
> On 26 October 2014 00:18, Hartmut Kaiser <har...@gm...> wrote: > At this point we assume, that polys[0] is a linear ring to be interpreted > as > a polygon exterior and polys[1:] are the corresponding interiors for > polys[0]. > > Here are our questions: > > Is this assumption correct? > Is there any detailed documentation describing the structure of the > returned > geometries? > Are the linear rings supposed to be correctly oriented (area > 0 for > exteriors and area < 0 for the interiors)? > > Hello Hartmut, > In brief, the answers are no, no and yes. > In more detail, assuming polys is not empty then it will contain one or > more polygon exteriors and zero or more interiors, and they can be in any > order. Here is a simple example where polys[0] is an interior and > polys[1] an exterior: > > x = [0, 0, 1, 1, 0.5] > y = [0, 1, 0, 1, 0.5] > z = [0.5, 0.5, 0.5, 0.5, 0] > triang = tri.Triangulation(x, y) > contour = plt.tricontourf(triang, z, levels=[0.2, 0.4]) > The returned geometries are purposefully not documented. They are an > 'implementation detail' and not considered part of the public interface. > and as such they could change at any time and hence should not be relied > upon. Of course you can choose to access them if you wish, as I do myself > sometimes, but we make no promises about what the order of the polygons > is, or that it won't change tomorrow. > In reality the order of the polygons is chosen to be something that is > easy for both the contour generation routines to create and for the > various backends to render. If you were to look at the output generated > by contourf, you will see that it is organised differently from that > produced by tricontourf and is more like you would like it to be, i.e. one > or more groups of an exterior polygon followed by zero or more > interiors. This is historical as the contourf code dates from before all > of the backends were able to render arbitrary groups of exterior and > interior polygons, and so the contourf code has to calculate the order for > the backends. When the tricontourf code was written the backends were all > able to calculate the exterior/interior containment themselves, so there > was no need for tricontourf to do it as well. Thanks Ian! Your detailed answer is much appreciated. As you might have already guessed, we have quite some problems creating clean geometries from the generated contour data. I have tried to put together one (reasonably) small test case illustrating at least one of our issues. I apologize for the lengthy code attached. The two attached figures demonstrate what we see. Matplotlib.png (generated by the library) does not really look ok. Also, the attached shape.png shows that there are a lot of geometries generated which are self-intersecting (highlighted in dark blue), and we already skip polygons with less than 3 points. BTW, there are many polygons stacked with the same geometries. Anything we can do about this? Thanks! Regards Hartmut --------------- http://boost-spirit.com http://stellar.cct.lsu.edu
On 1 November 2014 18:20, Hartmut Kaiser <har...@gm...> wrote: > Thanks Ian! Your detailed answer is much appreciated. > > As you might have already guessed, we have quite some problems creating > clean geometries from the generated contour data. I have tried to put > together one (reasonably) small test case illustrating at least one of our > issues. I apologize for the lengthy code attached. > > The two attached figures demonstrate what we see. Matplotlib.png > (generated by the library) does not really look ok. Also, the attached > shape.png shows that there are a lot of geometries generated which are > self-intersecting (highlighted in dark blue), and we already skip polygons > with less than 3 points. BTW, there are many polygons stacked with the same > geometries. > > Anything we can do about this? > > Thanks! > Regards Hartmut > Hartmut, You are using masked arrays where you shouldn't, again. The documentation for tricontour states that it expects z to be an array, it doesn't say masked array. If you pass it a masked array, it will ignore the mask. Hence you have a number of triangles that include a vertex with a z-value of -99999; when contoured these are going to give you lots of thin polygons that you don't want. You need to stop using masked arrays where they are not expected. Your triangulation should only contain triangles for which you have valid data at all three vertices. So either remove invalid triangles from your 'ele' array before creating the triangulation, or set a mask on the triangulation once it has been created, e.g. point_mask_indices = numpy.where(z.mask) tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1, 3), axis=1) triang.set_mask(tri_mask) Ian
Ian, > You are using masked arrays where you shouldn't, again. The documentation > for tricontour states that it expects z to be an array, it doesn't say > masked array. If you pass it a masked array, it will ignore the > mask. Hence you have a number of triangles that include a vertex with a > z-value of -99999; when contoured these are going to give you lots of thin > polygons that you don't want. > You need to stop using masked arrays where they are not expected. Your > triangulation should only contain triangles for which you have valid data > at all three vertices. So either remove invalid triangles from your 'ele' > array before creating the triangulation, or set a mask on the > triangulation once it has been created, e.g. > > point_mask_indices = numpy.where(z.mask) > tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1, > 3), axis=1) > triang.set_mask(tri_mask) Thanks very much for this explanation! With your code everything fell into place and we see the geometries we expect to see. Great library! Regards Hartmut --------------- http://boost-spirit.com http://stellar.cct.lsu.edu
Would it make sense to at least emit a warning when a mask is encountered. There are very few places in matplotlib where masked arrays are not allowed (I think histograms is the other spot, but I can't remember for sure). On Sun, Nov 2, 2014 at 11:10 AM, Ian Thomas <ian...@gm...> wrote: > On 1 November 2014 18:20, Hartmut Kaiser <har...@gm...> wrote: > >> Thanks Ian! Your detailed answer is much appreciated. >> >> As you might have already guessed, we have quite some problems creating >> clean geometries from the generated contour data. I have tried to put >> together one (reasonably) small test case illustrating at least one of our >> issues. I apologize for the lengthy code attached. >> >> The two attached figures demonstrate what we see. Matplotlib.png >> (generated by the library) does not really look ok. Also, the attached >> shape.png shows that there are a lot of geometries generated which are >> self-intersecting (highlighted in dark blue), and we already skip polygons >> with less than 3 points. BTW, there are many polygons stacked with the same >> geometries. >> >> Anything we can do about this? >> >> Thanks! >> Regards Hartmut >> > > Hartmut, > > You are using masked arrays where you shouldn't, again. The documentation > for tricontour states that it expects z to be an array, it doesn't say > masked array. If you pass it a masked array, it will ignore the mask. > Hence you have a number of triangles that include a vertex with a z-value > of -99999; when contoured these are going to give you lots of thin polygons > that you don't want. > > You need to stop using masked arrays where they are not expected. Your > triangulation should only contain triangles for which you have valid data > at all three vertices. So either remove invalid triangles from your 'ele' > array before creating the triangulation, or set a mask on the triangulation > once it has been created, e.g. > > point_mask_indices = numpy.where(z.mask) > tri_mask = numpy.any(numpy.in1d(ele, point_mask_indices).reshape(-1, > 3), axis=1) > triang.set_mask(tri_mask) > > Ian > > > ------------------------------------------------------------------------------ > > _______________________________________________ > Matplotlib-users mailing list > Mat...@li... > https://lists.sourceforge.net/lists/listinfo/matplotlib-users > >
On 2 November 2014 16:30, Benjamin Root <ben...@ou...> wrote: > Would it make sense to at least emit a warning when a mask is encountered. > There are very few places in matplotlib where masked arrays are not allowed > (I think histograms is the other spot, but I can't remember for sure). > Ben, That's certainly an option, but I'm happy with leaving it as it is with the onus on users to read the documentation. Ian