Polygon intersection area

Added by Nils about 1 month ago

Dear GMT users,

I have been trying to calculate the intersecting area of two polygons using GMT for a couple of days now, but I haven't been successful. The most promising result was with the gmt spatial -Si option. I know it says that it hasn't been implemented yet but it is doing something. Attached you'll find the two ASCII tables with the points of the polygons.

This is my code (simplified):

gmt gmtset PROJ_ELLIPSOID Sphere
gmt pscoast -R137:20/138:30/35:35/36:15 -JM137.476939/35.891485/20c -Ba20m/a20m -Df -G220/220/220 -W -Na -Ia/blue -A0/0/4 -Lf4/48.5/10/500+lkm -K > output.ps
# plot volcano
echo 137.476939 35.891485 '1.0' | gmt psxy -R137:20/138:30/35:35/36:15 -JM137.476939/35.891485/20c -St0.1i -Gblack -Wthinnest -O -K >> output.ps
# plot mapped deposit
gmt psxy Mapped_Deposit.txt -R137:20/138:30/35:35/36:15 -JM137.476939/35.891485/20c -L -W0.5,black,solid -O -K >> output.ps
# plot simulated deposit
gmt psxy Simulated_Deposit.txt -R137:20/138:30/35:35/36:15 -JM137.476939/35.891485/20c -L -W0.5,red,dashed -O -K >> output.ps
# find the intersecting polygon
gmt spatial Mapped_Deposit.txt Simulated_Deposit.txt -Si > Intersection.txt
# plot intersecting area
gmt psxy Intersection.txt -R137:20/138:30/35:35/36:15 -JM137.476939/35.891485/20c -L -W0.25,green,solid -O >> output.ps

I am also attaching the output.ps file. As you can see, there seem to be problems what the intersection points. Instead of following the path of the inner polygon it just jumps to the next intersection point. Does anybody have an idea on how to fix this, or maybe knows a different solution to my problem. In the end I just want to know the area of the intersecting polygon.

I am using GMT Version 5.4.2 (r18461).


Simulated_Deposit.txt Magnifier (54 KB)

Mapped_Deposit.txt Magnifier (78.3 KB)

output.ps (33.1 KB)

Replies (3)

RE: Polygon intersection area - Added by Peter about 1 month ago

Hi Nils

A crude method would be to use grdmask to produce two gridfiles with value 1 inside each polygon and 0 outside, then add the two gidfiles using grdmath and set all values <2 to 0 and then sum all cells, divide by two and multiply by the cell area (should be doable in a single call to grdmath).


RE: Polygon intersection area - Added by Nils about 1 month ago

Thanks Peter!
This certainly worked. I hope they soon implement the gmt spatial -Si option properly. It would be a lot faster.


RE: Polygon intersection area - Added by Joaquim about 1 month ago

Polygon algebra is not trivial to implement. At a time we talked about using the GEOS lib but it was never done.

You can do it in Mirone. The intersect area is 1057.258 km^2 (file attached)