We can do this fairly easily by converting our geographic regions into Region
s 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 FilledCurve
s which may be harder to turn into Region
s.
RegionIntersection
after makingRegions
from the polygons. $\endgroup$GeoDisk[geoPosition,Quantity[300,"Miles"]]
? $\endgroup$