More Advanced Data Cube Usage

Introduction to MODIS products and Snow Map

In this notebook you will learn more about the MODIS products in the Data Cube and how to use the Snow Cover and indices like 10-days NDSI.

This notebook will require a number of external Python modules to run, which will be imported at the top of the notebook as his convention.

As we have been showed in other notebooks, we need to initialise a Data Cube in order to have access to the Kyrgyzstan Data Cube (KDC), and this can be done as follows:

Now let's see which MODIS products are in the Data Cube using the command below (already showed in a previous training session).

As you can see there are many products in the Data Cube, let's see the details of one of the MODIS product. For example with the line below you can see which kind of data are in the product modis_nbar.

The results will tell you the data inside the product, the type of data, the nodata value, and the name with the aliases. Which means that you can load the data using either the measurement name or any of the aliases.

While with the following comand, you can check the description of a product. In this case 34 is the index of the product MOD10A1 found above from dc.list_products()

For MODIS the data have been downloaded from NASA, therefore, for more details, you can check the documentation online.

Now let's find the Snow Mask for the data in the Data Cube!

Once we have selected a Dataset, we can find how many cloud pixels are in our image. Rembember that the value for the cloudy pixel that we found from was 250

You can find the percentage of cloudy pixel in our raster. The .size function will give you the total number of pixels in our snow_cover raster

This means that 45% of our image is covered by clouds

Now we are computing how many snow pixels and nodata pixels are in our raster

Plotting our Snow Mask

We can load the shapefile of the raions using geopandas (see training session about geopandas)

The following command are used to plot the raster data (ds.snow_cover) within the geopandas dataframe raion_gdf

If you are interested in knowing only where the snow is, it is possible to use the snow_mask computed before. We can assign a value of 1 to all the pixels with snow, and 0 to the remaining ones.

If we want to do some statistical analysis with the raster data combined with our shapefile, we can use the module rasterstats. Click here for the documentation.

This module will allow us to operate zonal statistics. For example:

We copy all the results into a pandas DataFrame:

Now we have the count of the pixels for different classes (numbers in the header) and the total number of pixels a specific Raion (count). Then we add the data that we want into the geopandas DataFrame:

We can now plot the results in a nice looking map!

Percentage of snow by Raion:

Number of pixels that contains snow by Raion

Further analyses can be done using the snow cover.

In the following section we will look into snow persistence and snow percentage by pixel.

Snow Persistence and Snow Percentage

Let's get all the MODIS data for December 2019

Now we need to apply a NDSI limit to find our snow pixels

I group all my data by day, so I have a snow mask (snow/nosnow) for each day.

(The first dimension makes sense because there are 31 days in December.)

Now we can compute the snow persistence. This will tell us the maximum number of consecutive days where we had snow in a pixels.

We can also compute the snow percentage, which will tell us the number of days classified as snow over the number of observations (in this case, days)

First we compute the not valid pixels, where the pixels could not be classified as either snow or no_snow. For example, when the pixel has been identified as cloud, night or missing data.

Now we can compute the snow percentage