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:
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):
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
Save some space by converting to binary format: gmt convert land_polygons_osm_planet.gmt -bo2f > land_polygons_osm_planet.bf2
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.
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)
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:
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.
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.
