# How to display an intersection of two geo regions?

GeoGraphics[{
{
GeoStyling["StreetMap"],
First[GeoNearest["Country",geoPosition]]["Polygon"],GeoDisk[geoPosition,Quantity[300,"Miles"]]
}
GeoMarker[geoPosition]
},GeoBackground->None]


Suppose I evaluate the above expression with

geoPosition = GeoPosition[{35.6762,139.6503}]


It plots the background with a union of country's polygon and disk's polygon but I want intersection of the two. How to do this?

• I have more-or-less done this in this answer. I simply used RegionIntersection after making Regions from the polygons. Commented Feb 22, 2020 at 19:42
• @CarlLange do you know how to get polygon structure for the GeoDisk[geoPosition,Quantity[300,"Miles"]]? Commented Feb 22, 2020 at 20:11
• Easiest thing to do is generate a number of points equidistant from the centre and create a polygon from that. I'll add an answer in a moment. Commented Feb 22, 2020 at 20:40

We can do this fairly easily by converting our geographic regions into Regions and using RegionIntersection.

First let's get our geographic regions.

centre = GeoPosition[{35.6762, 139.6503}]

p1 = First[GeoNearest["Country", centre]]["Polygon"]


Now, it's not easy to use GeoDisk directly for our calculation, so we'll regenerate this by creating a Polygon with a list of points equidistant from the centre.

radius = Quantity[300, "Miles"]

p2 = Polygon[
Table[GeoDestination[centre, GeoDisplacement[{radius, b}]], {b, 0, 360,
1}]]


Now we will double-check that this looks right:

GeoGraphics[{p1, p2}]


Now we compute the intersection by turning the polygons into regions and taking RegionIntersection of that.

r1 = Region[p1 /. GeoPosition[x___] -> x]

r2 = Region[p2 /. GeoPosition[x___] -> x]

int = RegionIntersection[r1, r2, PerformanceGoal -> "Speed"]


Now we can turn this back into a geo polygon:

geoint = int /. Region[Polygon[x_]] -> Polygon[GeoPosition[x]]


And there we are:

GeoGraphics[{geoint}]


Now we can do things like GeoArea easily with GeoArea[geoint] (in this case, 178393km^2!)

Note that Sato island is not included - you may need to use GeoVariant[..., "AllAreas"] to get every extended part of a country. It may return FilledCurves which may be harder to turn into Regions.

• One more question, doesn't RegionIntersection also return a polygon object if the inputs are polygon? What is the purpose of MeshPrimitives[BoundaryDiscretizeRegion[...], 2]? Commented Feb 22, 2020 at 20:58
• Indeed, thank you, I have edited my answer. That was a holdover from some previous code I used as the base of this answer. Commented Feb 22, 2020 at 21:03
• Sorry I am wrong RegionIntersection returns MeshRegion in the case of polygon inputs so I think MeshPrimitives[..., 2] is necessary to get 2D-polygons. Commented Feb 22, 2020 at 21:59