Add a shapefiles tutorial

GenericMappingTools · Sep 2, 2023 · 7381a51 · 7381a51
1 parent 30b1f3a
commit 7381a51
Show file tree

Hide file tree

Showing 13 changed files with 187 additions and 11 deletions.
diff --git a/_assets/common_opts/opt_mfc.md b/_assets/common_opts/opt_mfc.md
@@ -1,4 +1,4 @@
-- **G** or **markerfacecolor** or **MarkerFaceColor** or **mc** or **fill**\
+- **G** or **markerfacecolor** or **MarkerFaceColor** or **markercolor** or **mc** or **fill**\
  Select color or pattern for filling of symbols [Default is no fill]. Note that plot will search for *fill*
  and *pen* settings in all the segment headers (when passing a GMTdaset or file of a multi-segment dataset)
  and let any values thus found over-ride the command line settings (but those must be provided in the terse GMT

diff --git a/documentation/common_features/color.md b/documentation/common_features/color.md
@@ -1,7 +1,7 @@
 # Setting color
 
-Color can be selected in several different ways. One of the is to create color maps with the *makecpt* and
-*grd2cpt* modules (see their own man pages). This is the method we use to colorize images, sets of points, etc.
+Color can be selected in several different ways. One of the is to create color maps with the \myreflink{makecpt} and
+\myreflink{grd2cpt} modules. This is the method we use to colorize images, sets of points, etc.
 The other option sets the color via keyword/value pairs and is appropriate to color fill polygons, individual
 symbols, etc and the one documented here.
 

diff --git a/documentation/modules/contour.md b/documentation/modules/contour.md
@@ -4,7 +4,7 @@
  contour(cmd0::String="", arg1=nothing; kwargs...)
 ```
 
-Contour table data by direct triangulation
+Contour plot from table data by direct triangulation
 
 
 Description

diff --git a/documentation/modules/grdcontour.md b/documentation/modules/grdcontour.md
@@ -4,9 +4,9 @@
 grdcontour(cmd0::String="", arg1=nothing; kwargs...)
 ```
 
-Make contour map using a grid.
+Make contour plot or map (using a projection) from a grid.
 
-Read a 2-D grid and produces a contour map by tracing each contour through the grid. Various options
+Read a 2-D grid and produces a contour plot by tracing each contour through the grid. Various options
 that affect the plotting are available. Alternatively, the *x, y, z* positions of the contour lines may be
 saved to one or more output files (or memory) and no plot is produced.
 

diff --git a/documentation/modules/scatter.md b/documentation/modules/scatter.md
@@ -27,6 +27,10 @@ Optional Arguments
 
 \textinput{common_opts/opt_mfc}
 
+- **zcolor** or **markerz** or **mz** : -- *zcolor=xx* **|** *zcolor=true*\
+ Take the vector `xx` (same size as number os points in data) and interpolate the current color scale to paint the
+ symbols based on that colr scale. The form `zcolor=true` is equivant to *zcolor=1:npoints*
+
 - **S** or *symbol* or *marker* or *Marker* or *shape* : -- Default is `circle` with a diameter of 7 points
  - *symbol=symbol* string\
  A full GMT compact string.

diff --git a/documentation/modules/triangulate.md b/documentation/modules/triangulate.md
@@ -18,11 +18,11 @@ starting at 0 for the first line) in the input file. As an option, you may choos
 multiple segment file that can be send to \myreflink{plot} to draw the triangulation network.
 If **region** and **inc** are set a grid will be calculated based on the surface defined by the
 planar triangles. The actual algorithm used in the triangulations is either that of Watson [1982]
-or Shewchuk [1996] [Default] (type **gmt("gmtget GMT_TRIANGULATE")** to see which method is
+or Shewchuk [1996] (Default. Type **gmt("gmtget GMT_TRIANGULATE")** to see which method is
 selected). This choice is made during the GMT (**not** `GMT.jl`) build step. Furthermore, if
 the Shewchuk algorithm is installed then you can also perform the calculation of Voronoi
 polygons and optionally grid your data via the natural nearest neighbor algorithm. **Note**:
-For geographic data with global or very large extent you should consider :doc:`sphtriangulate`
+For geographic data with global or very large extent you should consider \myreflink{sphtriangulate}
 instead since **triangulate** is a Cartesian or small-geographic area operator and is unaware
 of periodic or polar boundary conditions.
 

diff --git a/documentation/modules/trisurf.md b/documentation/modules/trisurf.md
@@ -6,7 +6,7 @@ trisurf(in, kw...)
 
 Plots the 3-D triangular surface defined by the points in a Mx3 matrix or a GMTdataset with data 
 x, y, z in the 3 first columns. The triangles are computed with a Delaunay triangulation done internaly.
-Since this is a \myreflink{plot3} _avatar_ all options in this function are those of the `plot3d` program.
+Since this is a \myreflink{plot3}  *avatar* all options in this function are those of the `plot3d` program.
 
 ---
 trisurf([x y z], kwargs...) plots the 3-D triangular surface defined by the points in vectors x, y, and z.

diff --git a/documentation/utilfuns.txt b/documentation/utilfuns.txt
@@ -24,6 +24,7 @@ image_alpha!
 image_cpt!
 imshow
 ind2rgb
+info
 inwhichpolygon
 isnodata
 kml2gmt
@@ -33,7 +34,8 @@ magic
 mat2ds
 mat2grid
 mat2img
-ODE2ds
+ODE2ds
+orbits
 pcolor
 plotgrid!
 plotyy

diff --git a/examples/contours.md b/examples/contours.md
@@ -1,3 +1,3 @@
-# Contours
+# Contour plots
 
 {{list_folder_with_images contours}}
diff --git a/tutorials.md b/tutorials.md
@@ -67,4 +67,15 @@
  ~~~</a>~~~
  @@
 
+ @@box
+ ~~~<a class="boxlink" href="vector_shp/vector_shp/">~~~
+ @@title Geospatial vector@@
+ @@box-content
+ ~~~
+ <img src="/tutorials/vector_shp/tilelogo.png">
+ ~~~
+ @@
+ ~~~</a>~~~
+ @@
+
 @@
diff --git a/tutorials/vector_shp.md b/tutorials/vector_shp.md
@@ -0,0 +1,3 @@
+# Reading, filtering and visualizing geospatial vector data
+
+{{list_folder_with_images vector_shp}}
diff --git a/tutorials/vector_shp/tilelogo.png b/tutorials/vector_shp/tilelogo.png
diff --git a/tutorials/vector_shp/vector_shp.md b/tutorials/vector_shp/vector_shp.md
@@ -0,0 +1,156 @@
+# Reading, filtering and visualizing geospatial vector data
+
+This tutorial is based on a [similar one for R](https://rpubs.com/geo2/vectordata) and shows how simpler and cleanly
+we can do this in Julia with GMT.jl. Vector data are composed of discrete geometric locations (*x, y* values) known
+as vertices that define the “shape” of the spatial object. There are three types of vector data: *Points, Lines* and
+*Polygons* (see more at, for example,
+[GIS in Python](https://www.earthdatascience.org/workshops/gis-open-source-python/intro-vector-data-python/)).
+Geospatial vector data is often stored in the shapefile format and that is the data format that we will use is this
+tutorial, but other formats could save been used without fundamental changes to the procedures shown below.
+
+## Reading vector data
+
+Read a shapefile representing Colombian departments which has been previously downloaded from DIVA-GIS,
+but not forgetting to load the package first.
+
+```julia
+using GMT
+deptos = gmtread(GMT.TESTSDIR * "COL_adm1.shp.zip");
+```
+
+We can see what the contents of the _deptos_ object are, but because there are many polygons we will restrict to show only the first 7 from the attribute table. For that we use the _info_ function that plots several types of information depending on the input data type. For GMTdatasets the _attribs_ keyword lets limit the number of table rows that are printed.
+
+
+```julia
+info(deptos, attribs=7)
+
+Attribute table (Dict{String, String})
+┌─────┬───────────┬──────┬─────────────┬───────────┬─────┬──────┬──────────────┬──────────┬───────────┐
+│ Row │ NAME_1 │ ID_1 │ ENGTYPE_1 │ NL_NAME_1 │ ISO │ ID_0 │ TYPE_1 │ NAME_0 │ VARNAME_1 │
+├─────┼───────────┼──────┼─────────────┼───────────┼─────┼──────┼──────────────┼──────────┼───────────┤
+│ 1 │ Amazonas │ 1 │ Commissiary │ │ COL │ 53 │ Comisaría │ Colombia │ │
+│ 2 │ Amazonas │ 1 │ Commissiary │ │ COL │ 53 │ Comisaría │ Colombia │ │
+│ 3 │ Antioquia │ 2 │ Department │ │ COL │ 53 │ Departamento │ Colombia │ │
+│ 4 │ Antioquia │ 2 │ Department │ │ COL │ 53 │ Departamento │ Colombia │ │
+│ 5 │ Antioquia │ 2 │ Department │ │ COL │ 53 │ Departamento │ Colombia │ │
+│ 6 │ Arauca │ 3 │ Intendancy │ │ COL │ 53 │ Intendencia │ Colombia │ │
+│ 7 │ Atlántico │ 4 │ Department │ │ COL │ 53 │ Departamento │ Colombia │ │
+└─────┴───────────┴──────┴─────────────┴───────────┴─────┴──────┴──────────────┴──────────┴───────────┘
+```
+
+It is useful to know the coordinate reference system of the vector data stored in the _deptos_ object. For that we will use again the _info_ function, but this time with the _crs_ keyword.
+
+```julia
+info(deptos, crs=true)
+
+PROJ: +proj=longlat +datum=WGS84 +no_defs
+WKT: GEOGCS["WGS 84",
+ DATUM["WGS_1984",
+ SPHEROID["WGS 84",6378137,298.257223563]],
+ PRIMEM["Greenwich",0],
+ UNIT["degree",0.0174532925199433]]
+```
+
+A fundamental starting point to understand data is to plot it. GMT can do it with great simplicity.
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+deptos = gmtread(GMT.TESTSDIR * "COL_adm1.shp.zip"); # Hide
+viz(deptos)
+```
+\end{examplefig}
+
+The somewhat large blank part on the upper part of the figure is due to small islands that are barely visible
+at this scale but also due to the automatic determination of the region boundaries. Also, since Colombia is
+close to equator the deformation implied by plotting geographical coordinates is not easily visible, but we
+can do better. We can tell GMT to guess a good projection for this and automatically apply it.
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+deptos = gmtread(GMT.TESTSDIR * "COL_adm1.shp.zip"); # Hide
+viz(deptos, proj=:guess)
+```
+\end{examplefig}
+
+## Filtering geospatial data based on attributes
+
+As we are interested only in one departament, we can filter the data. Review the following code and change
+it to match your department.
+
+```julia
+antioquia = filter(deptos, NAME_1 = "Antioquia");
+```
+
+Let’s plot the new object.
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+deptos = gmtread(GMT.TESTSDIR * "COL_adm1.shp.zip"); # Hide
+antioquia = filter(deptos, NAME_1 = "Antioquia"); # Hide
+viz(antioquia, proj=:guess)
+```
+\end{examplefig}
+
+We can repeat the previous steps to load the Colombian municipalities and filter the Antioquian ones.
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+munic = gmtread(GMT.TESTSDIR * "COL_ADM2.shp.zip");
+mun_antioquia = filter(munic, NAME_1=:Antioquia); # Symbols are as good as strings for attribute values
+viz(mun_antioquia, proj=:guess)
+```
+\end{examplefig}
+
+Now we are going to compute the centroids of the Antioquia municipes and use them as placers
+to contents of the *ID_2* attribute (that are numbers, though confess that I don't know what they represent).
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+munic = gmtread(GMT.TESTSDIR * "COL_ADM2.shp.zip"); # Hide
+mun_antioquia = filter(munic, NAME_1=:Antioquia); # Hide
+antioquia_points = centroid(mun_antioquia);
+t = info(mun_antioquia, att="ID_2"); # Get the values of the ID_2 attribute
+antioquia_points.text = t; # Add the tex column to the centroids object
+viz(mun_antioquia, proj=:guess, fill=:beige, lw=0, text=(data=antioquia_points, font=4))
+```
+\end{examplefig}
+
+And as an example, we can use the numeric information contained in the *ID_2* attribute (which comes out as string
+but we can easily convert it to numbers) and paint the municipe polygons according to those numbers. That is,
+make a choropleth map.
+
+\begin{examplefig}{}
+```julia
+using GMT # Hide
+munic = gmtread(GMT.TESTSDIR * "COL_ADM2.shp.zip"); # Hide
+mun_antioquia = filter(munic, NAME_1=:Antioquia); # Hide
+antioquia_points = centroid(mun_antioquia); # Hide
+t = info(mun_antioquia, att="ID_2"); # Hide
+# Convert t to numeric. Needed for creating a color map and make the choropleth style plot.
+tn = parse.(Int,t);
+
+# Create a color map
+C = makecpt(range=(minimum(tn),maximum(tn)), C=:bamako);
+
+# Vizualize 
+viz(mun_antioquia, proj=:guess, levels=tn, cmap=C, lw=0, title="Another Map of Antioquia", text=(data=antioquia_points, font=5), colorbar=true)
+```
+\end{examplefig}
+
+The examples shown here use vector data in the shapefile format, but the same could have been
+obtained using for example the geopackage format. In fact, as indicated in the SAGA-GIS site from
+where the shapefiles were downloaded, the original data originates from the [GADM](https://gadm.org)
+project where data was in _gpkg_ format.
+We can access the data from that site using the _gadm_ function (see examples in \myreflink{GADM}).
+For example, the Antioquia municipes can be downloaded with:
+
+```julia
+antioquia = gadm("COL", "Antioquia", children=true);
+```
+
+But that object does not contain the *ID_2* attribute so we couldn't in fact make the same choropleth map.