Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to use RegionIntersection in a couple of Polygons. The polygons are the following:

p1=Polygon[{{-102.251, 21.8628}, {-102.252, 21.8628}, {-102.252, 21.8631}, {-102.251, 21.8631}, {-102.251, 21.8628}}]

and

p2=Polygon[{{-102.253, 21.8671}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 21.8681}, {-102.252, 21.868}, {-102.252, 21.8678}, {-102.252, 21.8673}, {-102.253, 21.8672}, {-102.253, 21.8671}}]

When I run r = RegionIntersection[p1, p2] I was expecting that the resulting object was a "Region", however it is something like RegionIntersection[Polygon{{...}},Polygon{{...}}].

When I try to compute the Area of the resulting region it gives me:

Area::reg: "RegionIntersection[Polygon[{{-102.253,21.8671},{-102.253,21.8682},{-102.253,21.8682},{-102.253,21.8682},<<3>>,{-102.252,21.8673},{-102.253,21.8672},{-102.253,21.8671}}],<<1>>] is not a correctly specified region.

And if I test for a region with RegionQ[r] I obtain `true

Also if I try to plot the region discretizing with DiscretizeRegion[r] I obtain:

DiscretizeRegion[RegionIntersection[Polygon[{{-102.251, 21.8628}, {-102.252, 21.8628}, {-102.252, 21.8631}, {-102.251, 21.8631}, {-102.251, 21.8628}}],Polygon[{{-102.253, 21.8671}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 21.8681}, {-102.252, 21.868}, {-102.252, 21.8678}, {-102.252, 21.8673}, {-102.253, 21.8672}, {-102.253, 21.8671}}]]]

So it seems that the resulting object from the Intersection is a region but it is not possible to compute anything on it.

Following the examples on the Documentation I notice that using circles it behaves correctly and it is possible to compute on the intersection, however if the intersection is nothing (i.e. the two regions don't intersect) the object is still a region. I was wondering if there is a way of knowing that I have a null region. Of curse my first guess was computing the Area on it but for my polygons it is apparently impossible.

Any thoughts on why regions behave like that? and on how to know if the Intersection region is a null region?

share|improve this question
1  
The resulting object is a Region. It just turns out that the two polygons do not intersect and somehow they did not account for situations like this. You can see that there's no intersection by doing e.g. Graphics[{Red, p1, Blue, p2}]. It really should return a null region. –  RunnyKine Jul 17 at 4:08
    
I thought the same but I have similar behavior with RegionUnion –  moaimx Jul 17 at 4:26
2  
This looks like a bug and I filed it as such. –  user21 Jul 17 at 6:08
    
DiscretizeRegion seems to behave consistently on empty regions, but Area not working here is probably a bug. With RegionIntersection[Disk[{0, 0}, 1], Disk[{5, 5}, 1]] // Area you get 0. –  jtbandes Jul 17 at 7:52

1 Answer 1

up vote 2 down vote accepted

I think the problem is in the ordering of points in p2. if we plot p2 with the nodes we will find the following graphic:

p2 = Polygon[{{-102.253, 21.8671}, {-102.253, 21.8682}, {-102.253, 
     21.8682}, {-102.253, 21.8682}, {-102.253, 21.8681}, {-102.252, 
     21.868}, {-102.252, 21.8678}, {-102.252, 21.8673}, {-102.253, 
     21.8672}, {-102.253, 21.8671}}];
Show[Graphics[{Gray, p2}], 
 ListPlot[Tooltip /@ {{-102.253, 21.8671}, {-102.253, 
     21.8682}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 
     21.8681}, {-102.252, 21.868}, {-102.252, 21.8678}, {-102.252, 
     21.8673}, {-102.253, 21.8672}, {-102.253, 21.8671}}, 
  PlotStyle -> Directive[Red, PointSize[0.03]]]]

enter image description here

it can be seen that not all nodes are vertices of the polygon. if we remover those extra nodes the polygon will be like this:

p2 = Polygon[{{-102.253, 21.8681}, {-102.252, 21.868}, {-102.252, 
     21.8673}, {-102.253, 21.8672}}];
Show[Graphics[{Gray, p2}], 
 ListPlot[Tooltip /@ {{-102.253, 21.8681}, {-102.252, 
     21.868}, {-102.252, 21.8673}, {-102.253, 21.8672}}, 
  PlotStyle -> Directive[Red, PointSize[0.03]]]]

enter image description here

the area can be found easily like this:

    r = RegionUnion[p1, p2];
    Area[r]
   (* 1.1*10^-6*)

if you want to include all points you many rearrange them like this:

p2 = Polygon[{{-102.253, 21.8671}, {-102.252, 21.8673}, {-102.252, 
     21.8678}, {-102.252, 21.868}, {-102.253, 21.8682}, {-102.253, 
     21.8682}, {-102.253, 21.8682}, {-102.253, 21.8681}, {-102.253, 
     21.8672}, {-102.253, 21.8671}}];
Show[Graphics[{Gray, p2}], 
 ListPlot[Tooltip /@ {{-102.253, 21.8671}, {-102.252, 
     21.8673}, {-102.252, 21.8678}, {-102.252, 21.868}, {-102.253, 
     21.8682}, {-102.253, 21.8682}, {-102.253, 21.8682}, {-102.253, 
     21.8681}, {-102.253, 21.8672}, {-102.253, 21.8671}}, 
  PlotStyle -> Directive[Red, PointSize[0.03]]]]

enter image description here

and now the area can be found:

r = RegionUnion[p1, p2];
Area[r]
(*1.2*10^-6*)
share|improve this answer
    
I still think it's a bug; if you type out a simple Polygon that has points in the wrong order, you still can compute the area. –  jtbandes Jul 17 at 8:14
    
I will try it on other misbehaving polygons. Thanks! –  moaimx Jul 17 at 20:32

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.