<jupyter_start><jupyter_text>Data Structures and Export We will also learn about `xarray` - a data structure frequently used in geospatial analysis in Python, including in the Open Data Cube (ODC).Some of the content related to xarray in this notebook is derived from [Digital Earth Africa](https://www.digitalearthafrica.org/) notebooks: [here](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/master/Beginners_guide/07_Intro_to_xarray.ipynb) and [here](https://github.com/digitalearthafrica/deafrica-training-workshop/blob/master/docs/session_2/04_load_data_exercise.ipynb)  BackgroundWe will be discussing `Xarray`, which is a python library which simplifies working with labelled multi-dimension arrays. `Xarray` introduces labels in the forms of dimensions, coordinates and attributes on top of raw `numpy` arrays, allowing for more intitutive and concise development. We will also provide an overview of `pandas` `DataFrame` objects and `geopandas` `GeoDataFrame` objects, which are both tabular data structures, with a `GeoDataFrame` being a `DataFrame` with added functionality for working with shapes and shapefiles.This [external notebook](https://rabernat.github.io/research_computing/xarray.html) introduces more uses of xarray and may help you advance your skills further.More information about `xarray` data structures and functions can be found [here](http://xarray.pydata.org/en/stable/). More information about Numpy can be found [here](https://numpy.org/doc/stable/user/whatisnumpy.html).  DescriptionTopics covered include:* numpy data structures* xarray data structures* How to access xarray dimensions and metadata* Using indexing to explore multi-dimensional xarray data* Built-in xarray functions such as sum, std and mean* Plotting xarray data* pandas and geopandas* Writing data to storage***  Introduction to NumPyThe underlying data structure for all of the other data structures we will discuss in this section is Numpy arrays, or `ndarrays`. These arrays are multidimensional arrays of elements of the same type (e.g. integers, floating point numbers). The Numpy library has an assortment of routines for fast operations on arrays, including mathematical, logical, shape manipulation, sorting, selecting, I/O, discrete Fourier transforms, basic linear algebra, basic statistical operations, random simulation and much more.For more information on what NumPy is, see [here](https://numpy.org/doc/stable/user/whatisnumpy.html).NumPy arrays support slicing like lists do:<jupyter_code>import numpy as np
mylist = [1,2,3,4]
print(f'mylist[::2]: {mylist[::2]}')

myarr = np.array([1,2,3,4])
print(f'myarr[::2]: {myarr[::2]}')<jupyter_output>mylist[::2]: [1, 3]
myarr[::2]: [1 3]<jupyter_text>For multidimensional NumPy arrays, dimension indicies are separated by commas in indexing:<jupyter_code>list_2D = [[1,2,3], 
           [4,5,6]]
print(f'element 0,0 of list_2D: {list_2D[0][0]}')

arr_2D = np.array([[1,2,3], 
                   [4,5,6]])
print(f'element 0,0 of arr_2D: {arr_2D[0][0]}')<jupyter_output>element 0,0 of list_2D: 1
element 0,0 of arr_2D: 1<jupyter_text>Multidimensional NumPy arrays can be sliced, though nested lists cannot.<jupyter_code>print(f'first row of arr_2D: {arr_2D[:1]}')
print(f'first 2 elements of the first row of arr_2D: {arr_2D[:1,:2]}')<jupyter_output>first row of arr_2D: [[1 2 3]]
first 2 elements of the first row of arr_2D: [[1 2]]<jupyter_text>The shape of the array is available through the `.shape()` method of arrays.<jupyter_code>print(f'arr_2D: \n{arr_2D}')
print(f'arr_2D shape: \n{arr_2D.shape}')<jupyter_output>arr_2D: 
[[1 2 3]
 [4 5 6]]
arr_2D shape: 
(2, 3)<jupyter_text>Introduction to XarrayThe ODC uses xarray as its core data model. To better understand what it is, lets first do a simple experiment using a combination of plain numpy arrays and Python dictionaries.Suppose we have a satellite image with three bands: `Red`, `NIR` and `SWIR`. These bands are represented as 2-dimensional numpy arrays and the latitude and longitude coordinates for each dimension are represented using 1-dimensional arrays. Finally, we also have some metadata that comes with this image.  The code below creates fake satellite data and structures the data as a `dictionary`.<jupyter_code>#create faek satellite data
red = np.random.rand(250,250)
nir = np.random.rand(250,250)
swir = np.random.rand(250,250)

#create some lats and lons
lats = np.linspace(-23.5, -26.0, num=red.shape[0], endpoint=False)
lons = np.linspace(110.0, 112.5, num=red.shape[1], endpoint=False)

#create metadata
title = "Image of the desert"
date = "2019-11-10"

#stack into a dictionary
image = {"red": red,
         "nir": nir,
         "swir": swir,
         "latitude": lats,
         "longitude": lons,
         "title": title,
         "date": date}<jupyter_output><empty_output><jupyter_text>All our data is conveniently packed in a dictionary. Now we can use this dictionary to work with the data it contains:<jupyter_code>#date of satellite image
print(image["date"])

#mean of red values
image["red"].mean()<jupyter_output>2019-11-10<jupyter_text>Still, to select data we have to use numpy indexes. For example:<jupyter_code>image["red"][0,0]<jupyter_output><empty_output><jupyter_text>Wouldn't it be convenient to be able to select data from the images using the coordinates of the pixels instead of their relative positions? This is exactly what `xarray` solves! Let's see how it works:To explore xarray, we have a file containing some surface reflectance data. The object that we get is an xarray `Dataset`, which in some ways is very similar to the dictionary that we created before, but with lots of convenient functionality available.<jupyter_code>import xarray as xr

ds = xr.open_dataset('../resources/day4/data_structures_and_export/example_netcdf.nc')
ds<jupyter_output><empty_output><jupyter_text>Xarray dataset structureThere are 4 main elements of an `xarray.Dataset` object:* **Dimensions**: The dimensions of the dataset. For Earth observation data, this is often `x` (can be longitude), `y` (can be latitude) and `time`, as seen here. The units for the `x` and `y` dimensions are pixels, while `time` is counted in the number of scenes (images). In the example above, we see there are 12 images.* **Coordinates**: A list of the coordinate values of each dimension. `spatial_ref` refers to the [coordinate reference system (CRS)](https://en.wikipedia.org/wiki/Spatial_reference_system).* **Data variables**: The data values. We see `red`, `green` and `blue` are loaded in the example above. These values represent the amount of surface reflectance associated with these parts of the electromagnetic spectrum. These values are unitless.* **Attributes**: Metadata about this dataset. The CRS is listed again.A `Dataset` can be seen as a dictionary structure packing up the data, dimensions and attributes. Data variables in a `Dataset` object are called `DataArrays` and they share dimensions with the higher level `Dataset`. The figure below provides an illustrative example: A data variable of a `Dataset` can be accessed by its name through indexing notation, similarly to Python dictionaries. It can also be accessed by using `.` notation, which is more convenient.<jupyter_code>ds["green"]

#or alternatively 

ds.green<jupyter_output><empty_output><jupyter_text>There are attributes on this `DataArray` that apply uniquely to it (the data variable `green`), so are not on its containing `Dataset`:* **units**: The units for the values. `1` indicates that these values are unitless.* **nodata**: The value that denotes that a point's value is unknown.Dimensions are also stored as numeric arrays that we can easily access.<jupyter_code>ds['time']

#or alternatively 

ds.time<jupyter_output><empty_output><jupyter_text>Metadata is referred to as attributes and is internally stored under `.attrs`, but the same convenient `.` notation applies to them.<jupyter_code>ds.attrs['crs']

#or alternatively 

ds.crs<jupyter_output><empty_output><jupyter_text>DataArrays store their data internally as multidimensional NumPy arrays. But these arrays contain dimensions or labels that make it easier to handle the data. To access the underlaying NumPy array of a `DataArray` we can use the `.values` notation.<jupyter_code>arr = ds.green.values

type(arr), arr.shape<jupyter_output><empty_output><jupyter_text>IndexingXarray offers two different ways of selecting data. This includes the `isel()` approach, where data can be selected based on its index (like numpy).<jupyter_code>print(ds.time.values)

ss = ds.green.isel(time=0)
ss<jupyter_output>['2018-01-03T08:31:05.000000000' '2018-01-08T08:34:01.000000000'
 '2018-01-13T08:30:41.000000000' '2018-01-18T08:30:42.000000000'
 '2018-01-23T08:33:58.000000000' '2018-01-28T08:30:20.000000000'
 '2018-02-07T08:30:53.000000000' '2018-02-12T08:31:43.000000000'
 '2018-02-17T08:23:09.000000000' '2018-02-17T08:35:40.000000000'
 '2018-02-22T08:34:52.000000000' '2018-02-27T08:31:36.000000000']<jupyter_text>Or the `sel()` approach, used for selecting data based on its dimension of label value.<jupyter_code>ss = ds.green.sel(time='2018-01-08')
ss<jupyter_output><empty_output><jupyter_text>Slicing data is also used to select a subset of data.<jupyter_code>ss.x.values[100]
ss = ds.green.sel(time='2018-01-08', x=slice(2378390,2380390))
ss<jupyter_output><empty_output><jupyter_text>Built-in xarray functionsXarray exposes lots of functions to easily transform and analyse `Datasets` and `DataArrays`. > AggregationsAggregate values like the mean, standard deviation, or sum of an xarray can be calculated with similarly named method (`mean()`, `std()` and `sum()`, respectively).<jupyter_code>print("Mean of green band:", ds.green.mean().values)
print("Standard deviation of green band:", ds.green.std().values)
print("Sum of green band:", ds.green.sum().values)<jupyter_output>Mean of green band: 4141.488778766468
Standard deviation of green band: 3775.5536474649584
Sum of green band: 14426445446<jupyter_text>> Renaming data variablesSometimes the name of a data variable is not what we want it to be. For example, many data sources name data for the electromagnetic bands by their order for the data source rather than a familiar name like `green`. For example, for Landsat 8, the green band is band 3 ([source](https://www.usgs.gov/media/images/landsat-8-band-designations)), so may be named something like `B03`.To change the name of a data variable of a `Dataset`, use the `rename()` method.<jupyter_code>ds.rename({'blue': 'B02', 'green': 'B03', 'red': 'B04'})<jupyter_output><empty_output><jupyter_text>> MaskingSometimes certain values need to be set to NaN to avoid affecting computations (notably `nodata` values).The `where()` method can be used to mask undesired values. The first argument is an expression denoting which values should **not** be modified. The second argument is a value or array of the same shape to replace the masked values with. The second argument is optional and defaults to `NaN` so that xarray operations will ignore the masked points.Notice that, in this example, the mean after filtering `0` values (the `nodata` value for this data) is higher than before filtering.<jupyter_code>print(f'mean without filtering with where(): {ds.mean()}')
print(f'mean after filtering with where(): {ds.where(ds!=0).mean()}')<jupyter_output>mean without filtering with where(): <xarray.Dataset>
Dimensions:      ()
Coordinates:
    spatial_ref  int32 6933
Data variables:
    red          float64 3.984e+03
    green        float64 4.141e+03
    blue         float64 4.251e+03
mean after filtering with where(): <xarray.Dataset>
Dimensions:      ()
Coordinates:
    spatial_ref  int32 6933
Data variables:
    red          float64 4.292e+03
    green        float64 4.462e+03
    blue         float64 4.58e+03<jupyter_text>> Checking which values are within a rangeValues can be checked to be in a range with comparison operators like `` and operators like `&` (logical and). For example:<jupyter_code>print(f'mean of values in the range (4000, 5000): \n{ds.where((4000<ds) & (ds<5000)).mean()}')<jupyter_output>mean of values in the range (4000, 5000): 
<xarray.Dataset>
Dimensions:      ()
Coordinates:
    spatial_ref  int32 6933
Data variables:
    red          float64 4.464e+03
    green        float64 4.512e+03
    blue         float64 4.549e+03<jupyter_text>For more information about xarray, including more useful methods, see the documentation [here](https://xarray.pydata.org/en/stable/user-guide/index.html).  Plotting xarray dataPlotting is also conveniently integrated in the xarray library. `DataArray` objects have a `plot()` method that uses the `matplotlib` library to plot the data:<jupyter_code>ds["green"].isel(time=0).plot()<jupyter_output><empty_output><jupyter_text>but we still can do things manually using NumPy and Matplotlib if we choose:<jupyter_code>import matplotlib.pyplot as plt

rgb = np.dstack((ds.red.isel(time=0).values, 
                 ds.green.isel(time=0).values, 
                 ds.blue.isel(time=0).values))
rgb = np.clip(rgb, 0, 2000) / 2000

plt.imshow(rgb);<jupyter_output><empty_output><jupyter_text>But compare the above image to the result of elegantly chaining operations within xarray:<jupyter_code>ds[['red', 'green', 'blue']].isel(time=0).to_array().plot.imshow(robust=True, figsize=(6, 6));<jupyter_output><empty_output><jupyter_text>pandas and geopandasXarray helps us work with raster data, but sometimes we need to work with vector data (such as shapefiles) or tabular data.As mentioned previously, pandas `DataFrame` and geopandas `GeoDataFrame` objects are tabular data structures, with a `GeoDataFrame` being a `DataFrame` with added functionality for working with shapes and shapefiles.<jupyter_code>import geopandas as gpd

vector_file = '../resources/day4/data_structures_and_export/multiple_polygons_shapefile/multiple_polygons.shp'

gdf = gpd.read_file(vector_file)
gdf<jupyter_output><empty_output><jupyter_text>We often call the output `gdf` unless we have a better name. Sometimes the variable the object is written to ends with `_gdf` to clarify that it is a GeoDataFrame.The data in a `GeoDataFrame` can be visualized with the `head()` method. In this case, the result is the same as `gdf`, but `head()` only shows the first `n` rows, where `n` is its only argument (e.g. `head(10)` shows the first 10 rows). In contrast, `gdf` will show only a few rows - the first, last, and some in between, with the rest being represented by ellipses dots (`...`).<jupyter_code>gdf.head()<jupyter_output><empty_output><jupyter_text>For a `GeoDataFrame`, each row is a shape, stored in the `geometry` column. The shapes, denoted by `POLYGON (())` objects, have the type `shapely.geometry.polygon.Polygon`.A single record can be viewed in these ways:<jupyter_code>gdf.geometry.iloc[0]
# or
gdf.geometry[0]<jupyter_output><empty_output><jupyter_text>For more information about pandas, which applies to geopandas, see [this introduction to pandas](https://pandas.pydata.org/docs/user_guide/10min.html).For more information about geopandas and its uses, see [this documentation page](https://geopandas.org/getting_started/introduction.html).For information on loading data for the areas associated with shapes in a GeoDataFrame of a shapefile, see [this notebook](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/master/Frequently_used_code/Analyse_multiple_polygons.ipynb). For information on rasterizing a shapefile, which can be done with the [`dea_tools.spatial.xr_rasterize()`](https://docs.digitalearthafrica.org/en/latest/sandbox/gen/deafrica_tools.spatial.htmldeafrica_tools.spatial.xr_rasterize) function, see the “Rasterising a shapefile” section of [this notebook](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/master/Frequently_used_code/Rasterise_vectorise.ipynb).  Outputting Data > xarray to NetCDF**To export an xarray object to a NetCDF**, we use the `to_netcdf()` method of the xarray object:<jupyter_code>ds.to_netcdf('output.nc')<jupyter_output><empty_output><jupyter_text>> xarray to TIFF**To export an xarray.DataArray object to a TIFF**, we use the function [`datacube.utils.cog.write_cog()`](https://datacube-core.readthedocs.io/en/stable/dev/api/generate/datacube.utils.cog.write_cog.html). For an xarray called `ds` and a filepath string, such as 'output.tif', it is used as follows:<jupyter_code>from datacube.utils.cog import write_cog

write_cog(ds.red, fname='output.tif', overwrite=True)<jupyter_output><empty_output><jupyter_text>> xarray to vector (shapefile)**To export an xarray.DataArray object to a shapefile**, we use the function `dea_tools.spatialtools.xr_vectorize()`. For an xarray called `ds`, it is used as follows:<jupyter_code>import sys
import os
# Allow importing of our utilities.
sys.path.append(os.environ.get('NOTEBOOK_ROOT'))

from utils.dea_tools.spatialtools import xr_vectorize

gdf = xr_vectorize(ds.red)
gdf<jupyter_output><empty_output><jupyter_text>The result will be a `geopandas.GeoDataFrame` with a shape for each region of values. To export the result to a shapefile, specify the path of the shapefile to create in the `export_shp` parameter - which is a string path.For more information about about exporting data to TIFFs, see [this training document](https://github.com/digitalearthafrica/deafrica-sandbox-notebooks/blob/master/Frequently_used_code/Exporting_GeoTIFFs.ipynb). > pandas to CSV**To export a pandas object (including geopandas) to a shapefile**, we use the `to_csv()` method of the pandas object:<jupyter_code>gdf.to_csv('output.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10171790235081374, 'fraction_numerical': 0.036787974683544306, 'mean_word_length': 5.511037527593819, 'pattern_counts': {'":': 7, '<': 113, '<?xml version=': 0, '>': 119, 'https://': 16, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 12, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '29', 'dir': 'jupyter-structured-clean-dedup', 'id': '106062', 'license': 'Apache-2.0', 'n_tokens_mistral': 5965, 'n_tokens_neox': 5426, 'n_words': 1923, 'path': 'notebooks/day4/Data_Structures_and_Export.ipynb', 'repo_name': 'jcrattz/odc_training_notebooks', 'size': '17696'}