Update to the latest version

01-geodata.qmd

@@ -98,9 +98,9 @@ The concept of geospatial domain is formalized in
   general set of tools.
 - According to our definition of geospatial data, "raster data"
   is simply a table with colors as variables (e.g., RGB values)
-  combined with a Cartesian grid of quadrangle geometries.
-  We illustrate this concept in @fig-raster by zooming in a
-  satellite image of the Lena delta:
+  combined with a grid of quadrangle geometries. We illustrate
+  this concept in @fig-raster by zooming in a satellite image
+  of the Lena delta:
 
   ::: {#fig-raster layout-ncol=2}
 
@@ -268,9 +268,9 @@ The choice of table representation is a function of the application.
 
 All available domain representations come from the
 [Meshes.jl](https://github.com/JuliaGeometry/Meshes.jl)
-module. Let's start illustrating the most common
-domains in GIS, which are collections of disconnected
-2D geometries:
+module.
+
+Let's start with a simple list of disconnected geometries:
 
 ```{julia}
 using GeoStats
@@ -469,6 +469,14 @@ Similarly, we can add units to vectors of values:
 [1.0, 2.0, 3.0]u"m"
 ```
 
+It is also possible to load units explicitly to avoid the "u" prefix:
+
+```julia
+using Unitful: m, ft
+
+1.0m + 2.0ft
+```
+
 :::
 
 The function combines any table with any domain into a geospatial data
@@ -536,9 +544,16 @@ using the [GeoIO.jl](https://github.com/JuliaEarth/GeoIO.jl) module:
 ```{julia}
 using GeoIO
 
-GeoIO.load("data/geotable.geojson")
+GeoIO.load("data/countries.geojson")
 ```
 
+::: {.callout-note}
+
+The "data" folder is stored on GitHub.
+Check the [Preface](preface.qmd) for download instructions.
+
+:::
+
 We will see more examples of "vector data" in the chapter
 [Interfacing with GIS](03-geoio.qmd), and will explain why
 file formats like

02-geoviz.qmd

@@ -61,9 +61,9 @@ img = georef((A=rand(10, 10), B=rand(10, 10)))
 
 ::: {.callout-note}
 
-The `georef` function will create a `CartesianGrid` starting at the
-origin whenever the domain is omitted. The size of the grid is taken
-as the size of the first array in the named tuple:
+The `georef` function creates a `CartesianGrid` starting at the origin
+whenever the domain is omitted. The size of the grid is taken as the
+size of the first array in the named tuple:
 
 ```{julia}
 img.geometry
@@ -118,53 +118,36 @@ viz(img.geometry, color = img.B, alpha = rand(length(img.B)))
 ```
 
 Other aesthetic options are available in the official documentation.
-To really see the benefits of the framework, let's load data from a
-GeoJSON file and visualize it:
+
+As another example, consider the visualization of data over a `GeometrySet`:
 
 ```{julia}
-using GeoIO
+gset = GeometrySet([
+  Triangle((12, 12), (15, 15), (12, 15)),
+  Quadrangle((5, 12), (10, 15), (10, 18), (5, 15))
+])
 
-gis = GeoIO.load("data/geotable.geojson")
+gis = georef((A=[0.1, 0.2], B=[0.3, 0.4]), gset)
 
-viz(gis.geometry, color = 1:4)
+viz(gis.geometry, color = gis.A)
 ```
 
-::: {.callout-note}
-
-The "data" folder is stored on GitHub.
-Check the [Preface](preface.qmd) for download instructions.
-
-:::
-
-As already mentioned, the `viz!` function can be used to add more
-geometries to an existing scene. We can create a scene with the
-geometries from the first geotable ("raster data"), and then add
-the geometries from the second geotable ("vector data"):
+We can create a scene with the geometries from the first geotable ("raster data"),
+and then add the geometries from the second geotable ("vector data"):
 
 ```{julia}
-viz(img.geometry, color = 1:100)
-viz!(gis.geometry, color = 1:4)
+viz(img.geometry)
+viz!(gis.geometry)
 
 # display current figure
 Mke.current_figure()
 ```
 
-Let's add an additional set of points:
+Let's add an additional set of points and line segments to conclude the example:
 
 ```{julia}
-pts = [Point(-20, -10), Point(-20, 0), Point(-40, 10)]
-
-viz!(pts, color = 1:3)
-
-Mke.current_figure()
-```
-
-And a set of line segments to conclude the example:
-
-```{julia}
-seg = [Segment((-40, -10), (0, 0)), Segment((-40, 0), (-20, 10))]
-
-viz!(seg, color = 1:2)
+viz!([Point(-20, -10), Point(-20, 0), Point(-40, 10)])
+viz!([Segment((-40, -10), (0, 0)), Segment((-40, 0), (-20, 10))])
 
 Mke.current_figure()
 ```
@@ -220,6 +203,20 @@ fig
 
 :::
 
+The `viz` and `viz!` functions are aware of
+[Coordinate Reference Systems](https://en.wikipedia.org/wiki/Spatial_reference_system),
+which is a quite unique feature of the framework:
+
+```{julia}
+using GeoIO
+
+world = GeoIO.load("data/countries.geojson")
+
+viz(world.geometry, color = 1:nrow(world))
+```
+
+We will explore this concept more carefully in the chapter [Map projections](06-projections.qmd).
+
 ## viewer
 
 As geospatial data scientists we are often interested in quick inspection
@@ -229,20 +226,14 @@ To address this issue, the framework provides a basic `viewer` that displays
 all variables stored in a geotable:
 
 ```{julia}
-geotable = georef((A=rand(1000), B=rand(1000)), rand(3, 1000))
+geotable = georef((A=rand(1000), B=rand(1000)), rand(Point, 1000))
 
 viewer(geotable)
 ```
 
 ::: {.callout-note}
 
-The `georef` function will create a `PointSet` whenever its second argument
-is a `Matrix` of coordinates. In this case, the points are represented in the
-columns of the matrix:
-
-```{julia}
-geotable.geometry
-```
+The `rand(Point, 1000)` function call creates `1000` random `Point`s in 3D space.
 
 :::
 

03-geoio.qmd

@@ -85,18 +85,18 @@ various issues with available standards and backend implementations.
 Most GIS file formats do **not** preserve topological information.
 This means that neighborhood information is lost as soon as geometries
 are saved to disk. To illustrate this issue, we consider a geotable over
-a `CartesianGrid`:
+a grid:
 
 ```{julia}
-geotable = georef((A=rand(10, 10), B=rand(10, 10)))
+using GeoIO
+
+geotable = GeoIO.load("data/earth.tif")
 ```
 
 If we save the geotable to a `.geojson` file on disk, and then load it back,
 we observe that the `CartesianGrid` gets replaced by a `GeometrySet`:
 
 ```{julia}
-using GeoIO
-
 fname = tempname() * ".geojson"
 
 GeoIO.save(fname, geotable)
@@ -151,5 +151,10 @@ Notice how the `isconvex` function is compiled away to the **constant** `1` (i.e
 called on the triangle. The code for a generic polygon is much more complicated and
 requires runtime checks that are too expensive to afford, especially in 3D.
 
+Another reason to not adhere to a generic interface is that we can store information
+in the geometry types themselves (e.g., coordinate reference system) that is relevant to
+to design advanced scientific visualization features illustrated in the previous chapter,
+and to dispatch specialized algorithms from geodesic geometry.
+
 Having cleared that up, we will now proceed to the last foundational chapter of the book,
 which covers the advanced geometric processing features of the framework.

04-geoproc.qmd

@@ -129,8 +129,8 @@ polygonal areas in GIS:
 
 ```{julia}
 outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
-hole1 = [(0.2, 0.2), (0.4, 0.2), (0.4, 0.4), (0.2, 0.4)]
-hole2 = [(0.6, 0.2), (0.8, 0.2), (0.8, 0.4), (0.6, 0.4)]
+hole1 = [(0.2, 0.2), (0.2, 0.4), (0.4, 0.4), (0.4, 0.2)]
+hole2 = [(0.6, 0.2), (0.6, 0.4), (0.8, 0.4), (0.8, 0.2)]
 poly  = PolyArea([outer, hole1, hole2])
 
 viz(poly)
@@ -141,6 +141,29 @@ boundary of the `PolyArea`, also known as the outer `Ring`. The other two
 lists represent the two internal boundaries, or inner rings. A single list
 of vertices can be used, in which case the `PolyArea` doesn't have holes.
 
+::: callout-note
+
+The topology of a 2D `Polygon` is defined by the `orientation` of its rings.
+Most geometric processing algorithms expect counter clockwise (`CCW`) orientation
+for the outer ring, and clockwise (`CW`) orientation for the inner rings:
+
+```{julia}
+r = rings(poly)
+
+orientation.(r)
+```
+
+:::
+
+::: {.callout-note}
+
+## Tip for all users
+
+The `Repair` transform covered in the next chapter can be used to repair the
+orientations of rings in polygons.
+
+:::
+
 Finally, let's take a look into the polygonal `Chain`: 
 
 ```{julia}
@@ -151,8 +174,6 @@ These are 1-dimensional polytopes connecting `Point`s in sequence. We've
 seen the `Ring`s in the `PolyArea` and `Ngon` geometries:
 
 ```{julia}
-r = rings(poly)
-
 viz(r)
 ```
 
@@ -262,8 +283,8 @@ A very useful predicate is `intersects` (with a "s" at the end):
 
 ```{julia}
 outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
-hole1 = [(0.2, 0.2), (0.4, 0.2), (0.4, 0.4), (0.2, 0.4)]
-hole2 = [(0.6, 0.2), (0.8, 0.2), (0.8, 0.4), (0.6, 0.4)]
+hole1 = [(0.2, 0.2), (0.2, 0.4), (0.4, 0.4), (0.4, 0.2)]
+hole2 = [(0.6, 0.2), (0.6, 0.4), (0.8, 0.4), (0.8, 0.2)]
 poly  = PolyArea([outer, hole1, hole2])
 ball1 = Ball((0.5, 0.5), 0.05)
 ball2 = Ball((0.3, 0.3), 0.05)
@@ -376,9 +397,21 @@ for discretization, simplification, refinement, convex hull, etc.
 Below we illustrate some of these algorithms, which will be useful in future examples:
 
 ```{julia}
-points = rand(Point{2}, 100)
+points = rand(Point, 100, crs=Cartesian2D)
+```
+
+::: callout-note
+
+The `crs` option specifies the Coordinate Reference System (CRS) of the points.
+It will be explored in the chapter [Map projections](06-projections.qmd).
+By default, the `rand` function generates geometries with `Cartesian3D` crs:
+
+```{julia}
+rand(Triangle, 3)
 ```
 
+:::
+
 ```{julia}
 viz(boundingbox(points))
 viz!(points, color = "black")

05-transforms.qmd

@@ -525,8 +525,8 @@ let's consider the following polygonal area with holes:
 
 ```{julia}
 outer = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]
-hole1 = [(0.2, 0.2), (0.4, 0.2), (0.4, 0.4), (0.2, 0.4)]
-hole2 = [(0.6, 0.2), (0.8, 0.2), (0.8, 0.4), (0.6, 0.4)]
+hole1 = [(0.2, 0.2), (0.2, 0.4), (0.4, 0.4), (0.4, 0.2)]
+hole2 = [(0.6, 0.2), (0.6, 0.4), (0.8, 0.4), (0.8, 0.2)]
 poly  = PolyArea([outer, hole1, hole2])
 
 viz(poly)
@@ -618,11 +618,11 @@ the assignments, or request a grid size to discretize the `boundingbox` of all g
 ```{julia}
 A = [1, 2, 3, 4, 5]
 B = [1.1, 2.2, 3.3, 4.4, 5.5]
-p1 = PolyArea((2, 0), (6, 2), (2, 2))
-p2 = PolyArea((0, 6), (3, 8), (0, 10))
-p3 = PolyArea((3, 6), (9, 6), (9, 9), (6, 9))
-p4 = PolyArea((7, 0), (10, 0), (10, 4), (7, 4))
-p5 = PolyArea((1, 3), (5, 3), (6, 6), (3, 8), (0, 6))
+p1 = Triangle((2, 0), (6, 2), (2, 2))
+p2 = Triangle((0, 6), (3, 8), (0, 10))
+p3 = Quadrangle((3, 6), (9, 6), (9, 9), (6, 9))
+p4 = Quadrangle((7, 0), (10, 0), (10, 4), (7, 4))
+p5 = Pentagon((1, 3), (5, 3), (6, 6), (3, 8), (0, 6))
 gt = georef((; A, B), [p1, p2, p3, p4, p5])
 ```
 

07-pipelines.qmd

@@ -140,7 +140,7 @@ b = [1.6, 3.4, -1.0, 4.0]
 c = [3.4, 2.0, 3.6, -1.0]
 table = (; a, b, c)
 
-geotable = georef(table, rand(2, 4))
+geotable = georef(table, [(0, 0), (1, 0), (1, 1), (0, 1)])
 ```
 
 ```{julia}

08-splitcombine.qmd

@@ -175,18 +175,6 @@ geotable. To mark a string as a column name, we need to use curly braces:
 
 ::: {.callout-note}
 
-Alternatively, we could have used the following syntax to update the `geometry` column:
-
-```julia
-gt.geometry = box.(gt.geometry)
-```
-
-Currently, this syntax does **not** work with other columns.
-
-:::
-
-::: {.callout-note}
-
 ## Tip for all users
 
 The use of symbols to represent column names is preferred in macros.