How to display an intersection of two geo regions?

Question

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?

Answer 1

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.