Message Boards Message Boards

How to find the border coordinates

POSTED BY: Mohammad Bahrami
10 Replies

enter image description here - Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

POSTED BY: EDITORIAL BOARD

Hello Mads,

yes, yes, you are of course right: I did see it too simple! Iran-Azerbaijan: What a mean type of border this is! For my rehabilitation I might come up with this:

fbIR = Flatten[CountryData[Entity["Country", "Iran"], "FullCoordinates"], 1];
fbAZ = Flatten[CountryData[Entity["Country", "Azerbaijan"], "FullCoordinates"], 1];
isctn = Intersection[fbIR, fbAZ];
index = Map[First, 
  SplitBy[MapIndexed[{#1, #1 - First[#2]} &, 
    Sort@Flatten[Position[fbIR, #] & /@ isctn]], Last], {2}];  (* partitioning the point index list into connected parts *)
DynamicGeoGraphics[{Thick, GeoPath[GeoPosition /@ fbIR[[#]]] & /@ index}, ImageSize -> Large, Frame -> True]

But now the code becomes somewhat clumsy, and this I wanted to avoid ...

Regarding your second remark: The point list the GeoPath is plotted with must of course be the same as the index calculation is done with.

Many thanks for looking into my code and best regards -- Henrik

EDIT:

I just knew that there must be a special function for this: FindClusters is doing the trick!

fbIR = Flatten[CountryData[Entity["Country", "Iran"], "FullCoordinates"], 1];
fbAZ = Flatten[CountryData[Entity["Country", "Azerbaijan"], "FullCoordinates"], 1];
isctn = Intersection[fbIR, fbAZ];
index = FindClusters@Sort[Flatten[Position[fbIR, #] & /@ isctn]];(*partitioning the point index list into connected parts*)
DynamicGeoGraphics[{Thick, GeoPath[GeoPosition /@ fbIR[[#]]] & /@ index}, ImageSize -> Large, Frame -> True]
POSTED BY: Henrik Schachner

Thanks Henrik for the updated code! As I explained in the post, the boundary data can be very complex, and I tried to address these subtleties in my code.

Look at the case of India:

CountryData[Entity["Country", "India"], 
  "FullCoordinates"] // Dimensions

Let us try your code on India and Pakistan, and find the indexes wrt to India

fbIN = Flatten[
   CountryData[Entity["Country", "India"], "FullCoordinates"], 1];
fbPAK = Flatten[
   CountryData[Entity["Country", "Pakistan"], "FullCoordinates"], 1];
isctn = Intersection[fbIN, fbPAK];
index = FindClusters@
  Sort[Flatten[
    Position[fbIN, #] & /@ 
     isctn]];(*partitioning the point index list into connected \
parts*)DynamicGeoGraphics[{Thick, 
  GeoPath[GeoPosition /@ fbIN[[#]]] & /@ index}, ImageSize -> Large, 
 Frame -> True]

Note the ziczac feature on the left corner of the India-Pakistan border (close to the sea).

Also, as I said before, the result of your code is very sensitive to the country wrt which you find the indexes. For example, let's find indexes wrt to Pakisan:

fbIN = Flatten[
   CountryData[Entity["Country", "India"], "FullCoordinates"], 1];
fbPAK = Flatten[
   CountryData[Entity["Country", "Pakistan"], "FullCoordinates"], 1];
isctn = Intersection[fbIN, fbPAK];
index = FindClusters@
  Sort[Flatten[
    Position[fbPAK, #] & /@ 
     isctn]];(*partitioning the point index list into connected \
parts*)DynamicGeoGraphics[{Thick, 
  GeoPath[GeoPosition /@ fbIN[[#]]] & /@ index}, ImageSize -> Large, 
 Frame -> True]
POSTED BY: Mohammad Bahrami

Hi Mads,

again - thank you for looking into my code! Your second piece of code (in your last anwer) should read:

fbIN = Flatten[CountryData[Entity["Country", "India"], "FullCoordinates"], 1];
fbPAK = Flatten[CountryData[Entity["Country", "Pakistan"], "FullCoordinates"], 1];
isctn = Intersection[fbIN, fbPAK];
index = FindClusters@Sort[Flatten[Position[fbPAK, #] & /@ isctn]];(*partitioning the point index list into connected parts*)
DynamicGeoGraphics[{Thick, GeoPath[GeoPosition /@ fbPAK[[#]]] & /@ index}, ImageSize -> Large, Frame -> True]

i.e. the point list the GeoPath is plotted with must be the same as for the index calculation. And this in fact gives an seemingly good result. But this alone is of course not sufficient! So I do agree: My attempt is not robust. Consequently for calculating borders I will use your code! Again - many thanks for posting!

Best regards -- Henrik

As a final note: The border points at the location where my method breaks down, look pretty pathologic, it is difficult to see a border at all:

borderIndia = 
  CountryData[Entity["Country", "India"], "FullCoordinates"];
borderPakistan = 
  CountryData[Entity["Country", "Pakistan"], "FullCoordinates"];
Legended[DynamicGeoGraphics[{FaceForm[Red], 
   Entity["Country", "India"]["Polygon"], FaceForm[Green], 
   Entity["Country", "Pakistan"]["Polygon"], PointSize[.01],
   Green, Point[GeoPosition[borderPakistan]], PointSize[.005], Red, 
   Point[GeoPosition[borderIndia]]}], 
 SwatchLegend[{Green, Red}, {"Pakistan", "India"}]]

enter image description here

POSTED BY: Henrik Schachner

Thanks Henrik for your valuable comments! Although my previous code was working fine, I decided to add some revisions to my border function (motivated by our discussions). I am not very happy with the final code since it is not that simple; however, it addresses every form of complexities at the borders (like the odd case of Pakistan-India border).

BTW, I am not still sure what is going on with India's coordinates close to the sea at the border of Pakistan. I wonder if it is a computational bug in the Wolfram dataset or something else. Maybe @Jose M. Martin-Garcia can help us?

POSTED BY: Mohammad Bahrami

Hi Mads, excellent - many thanks! Best regards -- Henrik

POSTED BY: Henrik Schachner

Hello Mads,

Interesting problem, thanks for sharing! Given that Flatten is not changing the order, I think there is a somewhat simpler approach possible (if I am not seeing things too simple):

fbIR = Flatten[CountryData[Entity["Country", "Iran"], "FullCoordinates"], 1];
fbIQ = Flatten[CountryData[Entity["Country", "Iraq"], "FullCoordinates"], 1];
isctn = Intersection[fbIR, fbIQ];
minmax = MinMax@Flatten[Position[fbIR, #] & /@ isctn];
DynamicGeoGraphics[{Thick, GeoPath[GeoPosition /@ fbIR[[Span @@ minmax]]]}, ImageSize -> Large,Frame -> True]

Using this very nice DynamicGeoGraphics one can see that the border points are still somewhat coarse.

Regards -- Henrik

POSTED BY: Henrik Schachner

Thanks Henrik! Neat! However, borders are sometimes disconnected parts and your approach won't work in those cases. For example:

fbIR = Flatten[
   CountryData[Entity["Country", "Iran"], "FullCoordinates"], 1];
fbAZ = Flatten[
   CountryData[Entity["Country", "Azerbaijan"], "FullCoordinates"], 1];
isctn = Intersection[fbIR, fbAZ];
minmax = MinMax@Flatten[Position[fbIR, #] & /@ isctn]
DynamicGeoGraphics[{Thick, 
  GeoPath[GeoPosition /@ fbIR[[Span @@ minmax]]]}, ImageSize -> Large,
  Frame -> True]

Additionally, an important question is which country should be used in computing minmax? In case of Iran-Azerbaijan, if one uses fbAZ, the result will be very off.

fbIR = Flatten[
   CountryData[Entity["Country", "Iran"], "FullCoordinates"], 1];
fbAZ = Flatten[
   CountryData[Entity["Country", "Azerbaijan"], "FullCoordinates"], 1];
isctn = Intersection[fbIR, fbAZ];
minmax = MinMax@Flatten[Position[fbAZ, #] & /@ isctn]
DynamicGeoGraphics[{Thick, 
  GeoPath[GeoPosition /@ fbIR[[Span @@ minmax]]]}, ImageSize -> Large,
  Frame -> True]
POSTED BY: Mohammad Bahrami

Dear Mads,

Thanks for posting. I really like this; it is useful for some projects I run with students. There are interesting problems about the fractal dimensions of countries and coastlines, or adjacent countries reporting different lengths of the same common border.

Thanks,

Marco

POSTED BY: Marco Thiel

Thanks Marco! That would be great if you share some of those problems :-)

POSTED BY: Mohammad Bahrami
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract