Adding code lines

02-Spaces.qmd

@@ -248,19 +248,61 @@ tradition of specifying points as $(\phi,\lambda)$ but throughout
 this book we will stick to longitude-latitude, $(\lambda,\phi)$.
 The point denoted in @fig-sphere  has $(\lambda,\phi)$ or ellipsoidal
 coordinates with angular values
+
+::: panel-tabset
+
+#### R
 ```{r echo=!knitr::is_latex_output()}
 #| code-fold: true
 #| collapse: false
 p <- st_as_sfc("POINT(60 47)", crs = 'OGC:CRS84')
 p[[1]]
 ```
+#### Python
+```{python}
+#| code-fold: true
+#| collapse: false
+import geopandas as gpd
+from shapely.geometry import Point
+
+point = Point(60, 47)
+crs = 'epsg:4326'
+
+gdf = gpd.GeoDataFrame(geometry=[point], crs=crs)
+p = gdf.geometry.values[0]
+print(p)
+
+```
+:::
 measured in degrees, and geocentric coordinates
+
+::: panel-tabset
+
+#### R
 ```{r echo=!knitr::is_latex_output()}
 #| code-fold: true
 #| collapse: false
 p <- st_as_sfc("POINT(60 47)", crs = 'OGC:CRS84') |> st_transform("+proj=geocent")
 p[[1]]
 ```
+#### Python
+```{python}
+#| code-fold: true
+#| collapse: false
+import geopandas as gpd
+from shapely.geometry import Point
+import pyproj
+
+point = Point(60, 47)
+crs = 'epsg:4326'  # assuming OGC:CRS84 is equivalent to EPSG:4326
+
+gdf = gpd.GeoDataFrame(geometry=[point], crs=crs)
+projected_crs = pyproj.CRS.from_proj4("+proj=geocent")
+gdf_projected = gdf.to_crs(projected_crs)
+p = gdf_projected.geometry.values[0]
+print(p)
+```
+:::
 measured in metres.
 
 \index{coordinates!units}

03-Geometries.qmd

@@ -384,12 +384,33 @@ and it represents essentially the empty set: when combining
 it vanishes.
 
 All geometry types have a special value representing the empty (typed) geometry, like
+
+::: panel-tabset
+
+#### R
 ```{r echo=!knitr::is_latex_output()}
 #| code-fold: true
 #| collapse: false
 st_point()
 st_linestring(matrix(1,1,3)[0,], dim = "XYM")
 ```
+#### Python
+```{python}
+#| code-fold: true
+#| collapse: false
+import geopandas as gpd
+from shapely.geometry import Point
+from shapely.geometry import LineString
+
+point = Point()
+linestring = LineString()
+gdf_point = gpd.GeoDataFrame(geometry=[point])
+gdf_line = gpd.GeoDataFrame(geometry=[linestring])
+print(gdf_point)
+print(gdf_line)
+```
+
+:::
 and so on, but they all point to the empty set, differing only in their
 dimension (@sec-de9im).
 

09-Large.qmd

@@ -58,13 +58,31 @@ with a WKT text string containing a geometry; only geometries from
 the target file that intersect with this geometry will be returned.
 An example is
 
+::: panel-tabset
+
+#### R
 ```{r}
 library(sf)
 file <- system.file("gpkg/nc.gpkg", package = "sf")
 c(xmin = -82,ymin = 36, xmax = -80, ymax = 37) |>
     st_bbox() |> st_as_sfc() |> st_as_text() -> bb
 read_sf(file, wkt_filter = bb) |> nrow() # out of 100
 ```
+#### Python
+```{python}
+import geopandas as gpd
+from shapely.geometry import box
+
+# Define the bounding box and read the data
+bbox = box(-82, 36, -80, 37)
+gdf = gpd.read_file("https://github.com/r-spatial/sf/raw/main/inst/gpkg/nc.gpkg")
+gdf_bbox = gdf[gdf.intersects(bbox)]
+# Print the number of rows in the data
+print(gdf_bbox.shape[0])
+
+```
+
+:::
 
 Here, `read_sf` is used as an alternative to `st_read` to suppress
 output.

10-Models.qmd

@@ -69,12 +69,39 @@ Objects of class `sf` need no special treatment, as they are
 `data.frame`s. To create maps of the resulting predictions, predicted
 values need to be added to the `sf` object, which can be done using the `nc`
 dataset loaded as in @sec-intro by:
+
+::: panel-tabset
+
+#### R
 ```{r sid, cache = FALSE}
 nc |> mutate(SID = SID74/BIR74, NWB = NWBIR74/BIR74) -> nc1
 lm(SID ~ NWB, nc1) |>
   predict(nc1, interval = "prediction") -> pr
 bind_cols(nc, pr) |> names()
 ```
+#### Python
+```{python}
+import geopandas as gpd
+import statsmodels.formula.api as smf
+
+# Read the 'nc' dataset from the GeoPackage file
+file = gpd.read_file("https://github.com/r-spatial/sf/raw/main/inst/gpkg/nc.gpkg")
+
+# Create the new columns and fit the linear model
+nc1 = file.assign(SID=file["SID74"]/file["BIR74"], NWB=file["NWBIR74"]/file["BIR74"])
+model = smf.ols(formula="SID ~ NWB", data=nc1).fit()
+
+# Compute the predictions and prediction intervals
+pr = model.get_prediction(nc1).summary_frame()
+
+# Combine the original data with the predictions
+result = pd.concat([file, pr], axis=1)
+
+# Print the column names
+print(result.columns)
+```
+
+:::
 where we see that 
 
 * `lm` estimates a linear model and works directly on an `sf` object

11-PointPattern.qmd

@@ -60,25 +60,68 @@ include
 
 Point patterns have an observation window. Consider the points generated randomly by 
 
+::: panel-tabset
+
+#### R
 ```{r}
 library(sf)
 n <- 30
 set.seed(13531) # remove this to create another random sequence
 xy <- data.frame(x = runif(n), y = runif(n)) |> 
     st_as_sf(coords = c("x", "y"))
 ```
+#### Python
+```{python}
+import geopandas as gpd
+from shapely.geometry import Point
+import random
+
+random.seed(13531)
+n = 30
+x = [random.uniform(0, 1) for i in range(n)]
+y = [random.uniform(0, 1) for i in range(n)]
+xy = gpd.GeoDataFrame(geometry=[Point(x, y) for x, y in zip(x, y)])
+xy.crs = "EPSG:4326"
+
+# Print the GeoDataFrame
+print(xy)
+
+```
+:::
 
 then these points are (by construction) uniformly distributed,
 or completely spatially random, over the domain $[0,1] \times [0,1]$. For
 a larger domain, they are not uniform, for the two square windows
 `w1` and `w2` created by 
+
+::: panel-tabset
+
+#### R
 ```{r}
 w1 <- st_bbox(c(xmin = 0, ymin = 0, xmax = 1, ymax = 1)) |> 
 		st_as_sfc() 
 w2 <- st_sfc(st_point(c(1, 0.5))) |> st_buffer(1.2)
 ```
+#### Python
+
+```{python}
+import geopandas as gpd
+from shapely.geometry import box, Point
+
+bbox = box(0, 0, 1, 1)
+w1 = gpd.GeoDataFrame(geometry=[bbox])
+point = Point(1, 0.5)
+buffer = point.buffer(1.2)
+w2 = gpd.GeoDataFrame(geometry=[buffer])
+print(w1, w2)
+```
+:::
+
 this is shown in @fig-pp1.
 
+::: panel-tabset
+
+#### R
 ```{r fig-pp1, echo=!knitr::is_latex_output()}
 #| fig.height: 3.3
 #| fig.cap: "Depending on the observation window (grey), the same point pattern can appear completely spatially random (left), or clustered (right)"
@@ -89,17 +132,58 @@ plot(xy, add = TRUE)
 plot(w2, axes = TRUE, col = 'grey')
 plot(xy, add = TRUE, cex = .5)
 ```
+#### Python
+
+```{python}
+#| fig.height: 3.3
+#| fig.cap: "Depending on the observation window (grey), the same point pattern can appear completely spatially random (left), or clustered (right)"
+#| code-fold: true
+import geopandas as gpd
+import matplotlib.pyplot as plt
+
+# Set the plot parameters
+fig, axs = plt.subplots(1, 2, figsize=(10, 5))
+axs[0].set_xlim([0, 1])
+axs[0].set_ylim([0, 1])
+axs[1].set_xlim([-0.2, 2.2])
+axs[1].set_ylim([-0.2, 1.2])
+
+# Plot the first map
+w1.plot(ax=axs[0], color='grey', alpha=0.5, edgecolor='black')
+xy.plot(ax=axs[0], markersize=50, color='red')
+
+# Plot the second map
+w2.plot(ax=axs[1], color='grey', alpha=0.5, edgecolor='black')
+xy.plot(ax=axs[1], markersize=25, color='red')
+
+# Show the plot
+plt.show()
+
+```
+:::
+
 
 \index{ppp}
 \index{sf!as.ppp}
 
 Point patterns in **spatstat** are objects of class `ppp` that contain
 points and an observation window (an object of class `owin`).
 We can create a `ppp` from points by
+
+::: panel-tabset
+#### R
 ```{r}
 library(spatstat) |> suppressPackageStartupMessages()
 as.ppp(xy)
 ```
+#### Python
+```{python}
+from pointpats import PointPattern
+coords = [(point.x, point.y) for point in xy.geometry]
+p1 = PointPattern(coords)
+p1.summary()
+```
+:::
 where we see that the bounding box of the points is used as observation
 window when no window is specified. If we add a polygonal geometry as the
 first feature of the dataset, then this is used as observation window: