simple raster write example

pygis/docs/c_new_vectors.md

@@ -32,7 +32,7 @@ html_meta:
 ----------------
 
 
-# Creating Spatial Vector Data
+# Spatial Points Lines Polygons in Python
 We often find ourselves in a situation where we need to generate new spatial data from scratch, or need to better understand how our data is constructed. This lesson will walk you through some of the most common forms of data generation. 
 ```{code-cell} ipython3
 # Import necessary modules first

pygis/docs/c_rasters.md

@@ -28,7 +28,7 @@ html_meta:
 ```
 ----------------
 
-# Raster Data 
+# Raster Data in Python
 A raster data model uses an array of cells, or pixels, to represent real-world objects. Raster datasets are commonly used for representing and managing imagery, surface temperatures, digital elevation models, and numerous other entities. A raster can be thought of as a special case of an area object where the area is divided into a regular grid of cells. But a regularly spaced array of marked points may be a better analogy since rasters are stored as an array of values where each cell is defined by a single coordinate pair inside of most GIS environments. Implicit in a raster data model is a value associated with each cell or pixel. This is in contrast to a vector model that may or may not have a value associated with the geometric primitive.
 
 In order to work with raster data we will be using `rasterio` and later `geowombat`. Behind the scenes a `numpy.ndarray` does all the heavy lifting. To understand how raster works it helps to construct one from scratch. 
@@ -59,3 +59,56 @@ plt.imshow(Z)
 plt.title("Temperature")
 plt.show()
 ``` 
+Note that *`Z` contains no data on its location*. Its just an array, the information stored in `x` and `y` aren't associated with it at all. This location data will often be stored in the header of file. In order to 'locate' the array on the map we will use affine transformations. 
+
+
+
+## Assigning spatial data to an array in python
+Ok we have an array of data and some coordinates for each cell, but how do we create a spatial dataset from it? In order to do this we need three components:
+
+1) An array of data and typically the xy coordinates
+2) A coordinate reference system which defines what coordinate space we are using (e.g. degrees or meters, where is the prime meridian etc)
+3) A transform defining the coordinate of the upper left hand corner and the cell resolution
+
+```{note}
+These topic is covered extensively in the next chapter. We will briefly introduce the topic here, but will go into the details later. For more info on transforms go [here](d_raster_crs_intro.md). For more info on coordinate references systems go [here](d_crs_what_is_it.md) or to understand CRS codes go [here](https://pygis.io/docs/d_understand_crs_codes.html).
+```
+
+Once you have those components you can write out a working spatial raster dataset in python in a few lines of code. We just need to provide the information listed above in a format that rasterio understands. 
+
+
+ ```{code-cell} ipython3
+from rasterio.transform import Affine
+import rasterio as rio
+
+res = (x[-1] - x[0]) / 240.0
+transform = Affine.translation(x[0] - res / 2, y[0] - res / 2) * Affine.scale(res, res)
+
+# open in 'write' mode, unpack profile info to dst
+with rio.open(
+    "../temp/new_raster.tif",
+    "w",
+    driver="GTiff",   #output file type
+    height=Z.shape[0], # shape of array
+    width=Z.shape[1],
+    dtype=Z.dtype,  # output datatype
+    crs="+proj=latlong",  # CRS
+    transform=transform,  #location and resolution of upper left cell
+) as dst:
+
+    # check for multiband
+    if len(Z.shape) == 3:
+        # write each band individually
+        for band in range(len(Z)):
+            # write data, band # (starting from 1)
+            dst.write(Z[band], band + 1)
+    # write single band
+    else:
+        dst.write(Z, 1)
+```
+
+```{note}
+Raster data is often 'multiband' (e.g. red, green, blue), so I provided code that works for both multiband and single band data. 
+
+If you are storing multiband data the dimensions should be stored as (band, y, x ).
+```