Trouble with open polygons on plot borders

I’m using the OpenStreetMap (OSM) coast lines as described in ex51 for a plot of Prince Edward Island in Canada.

Everything works fine when using the OSM global data set coming in at a whopping ~568MB in binary form. Executing the script took ~14s on my machine:

Cutting out the area of interest to speed things up reduces the binary file size to ~0.5MB and the execution time of the script to < 1s (similar to GSHHG). But I get in trouble at the lower left corner of the plot:

Nova Scotia and New Brunswick are heavily messed up due to – what I believe – are open polygons at the plot border.

My question: Is there an easy way to prevent this while maintaining small file sizes and quick execution times?

Some stuff to play with:

Generate the second plot with the messed up lower left corner: (239.3 KB)

gmt begin map_Prince_Edward_Island_OSM_cutout png
  gmt plot -R-64.6/45.8/-61.7/47.2+r -JL-63.15/46.5/46.2/46.8/20c -Bag \
    -B+t"Prince Edward Island - Canada"+s"OpenStreetMap cutout"+glightcyan1 \
    land_polygons_osm_pei.bf2 -bi2f -Wthinnest,cadetblue1 -Ggrey
gmt end show

Get the global OSM coast lines and make the correct plot (GDAL required):

  1. Download the OSM coastlines shape file (Caution! ZIP file size 808.4 MB)

  2. Convert to something useful (Caution! This will generate a file of ~1.6 GB):
    ogr2ogr -f OGR_GMT land_polygons_osm_planet.gmt land_polygons.shp

  3. Save some space by converting to binary format:
    gmt convert land_polygons_osm_planet.gmt -bo2f > land_polygons_osm_planet.bf2

  4. Make the same plot as above but with Nova Scotia and New Brunswick intact:


gmt begin map_Prince_Edward_Island_OSM_planet png
  gmt plot -R-64.6/45.8/-61.7/47.2+r -JL-63.15/46.5/46.2/46.8/20c -Bag \
    -B+t"Prince Edward Island - Canada"+s"OpenStreetMap planet"+glightcyan1 \
    land_polygons_osm_planet.bf2 -bi2f -Wthinnest,cadetblue1 -Ggrey
gmt end show

First, the way you want to select your data is clipping to the area of interest rather than cutting. Clipping is supposed to produce closed polygons, while cutting only selects those points that are inside your region of interest. I don’t know how to clipping polygon data to another polygon using GMT commands, not even sure this is possible.

I have done polygon data clipping myself using ogr2ogr with -clipsrc option. However ogr2ogr is not supposed to recognize GMTs binary formats (those specified using -bo2f and -bi2f) as far as I can understand. It can however work with the original OSM coastline shapefile.

Second, a minor point. You should be careful when choosing an appropriate region of interest. Like you choose a rectangular area using geographical coords, while your plotting script uses a combination of projection and region that is not rectangular with regard to geographical coords. But you will surely see if your ROI is smaller than the mapping region.

1 Like

I can make it work with GMT.jl but it’s slow
I think it should be also possible (and wiser) to do what Mikhail advised but did not try it.

# This is fast
Dosm = gmtread("C:\\v\\land-polygons-complete\\land_polygons.shp");

# This very slow. Probably because it is not using any index file
D2 = intersection(Dosm, mat2ds([-64.5 45.8; -64.5 47.08; -61.8 47.08; -61.8 45.8; -64.5 45.8], geom=3))

plot(D2, G=:gray,  show=1)

Hello @mkononets and @Joaquim, a sincere thank you for your answers.

I tried Joaquim’s GMT.jl variant as I had Julia already running. It looks way better than my plot! A bit slow maybe. Unfortunately it eats Sheopdy Bay and Cumberland Basin – the waters in the lower left corner of the plot bordering the map frame.

Using Mikhail’s suggestion with ogr2ogr -clipsrc worked beautifully:

The commands I used:

ogr2ogr -f OGR_GMT land_polygons_osm_pei_clipped.gmt land_polygons.shp -clipsrc -65 45 -61 48

gmt convert land_polygons_osm_pei_clipped.gmt -bo2f > land_polygons_osm_pei_clipped.bf2

gmt begin map_Prince_Edward_Island png
  gmt plot -R-64.6/45.8/-61.7/47.2+r -JL-63.15/46.5/46.2/46.8/20c -Bag \
    -B+t"Prince Edward Island - Canada"+s"ogr2ogr -clipsrc"+glightcyan1 \
    land_polygons_osm_pei_clipped.bf2 -bi2f -Wthinnest,cadetblue1 -Ggrey
gmt end show

The ogr2ogr command took the longest time to execute with ~19s. The call to gmt convert is not needed performance wise. The plot is generated in <1s for both the binary file and the ASCII text file. For storing the text file is roughly three times bigger than the binary file.

Thank you also for bringing the difference between cutting and clipping a polygon to my attention. It does make a crucial difference when working with closed polygons!

Additionally you are totally right about my reckless way of selecting my region of interest – I was too lazy making it a bit bigger and got away with it this time. I even wrote a feature request some years ago to support my laziness. :wink:

If you use the same limits in the ogr2ogr as I used in the intersection you’ll get a plot exactly like mine.

For some annoying reason calling ogr2ogr from Julia is not clipping to the limits but writes the entire file contents. And take ages if I ask it to return data to Julia instead of saving to disk plus some other oddities.

Just a little bit off-topic: That GMT.jl package is awesome @Joaquim !

1 Like