TopoSUB clustering with land-cover types for areas of irregular shapes

[achillejouberton]

Dear TopoPyScale community,

My name is Achille Jouberton, a PhD student within Francesca Pellicciotti’s group at WSL (Zürich, Switzerland).
My research focuses on land-surface modelling at glacierized catchments in High Mountain Asia, focusing on snowfall dynamics and glacier mass balance. I am using a land-surface model at 100m spatial resolution over catchments which are 100-300 km^2 in size, and since a couple of months I have attempted to use the clustering approach implemented in TopoPyScale to perform much faster fully distributed simulations (yet neglecting lateral exchanges). The results are quite promising, I show an example of such application at Parlung No.4 glacier (Fig. 1), the results between the clustered approach (T&C-175-Clusters) and the fully distributed run (T&C) are very similar for snow cover fraction (Fig. 2) and distributed glacier surface mass balance (Fig.3). The computational time went down from ~35 hours to ~30min by using the cluster approach.
Using TopoPyScale, I faced the following challenges and came up with the following temporary solutions whenever necessary:

1. TopoPyScale seems to only accept a complete – rectangular – DEM with no missing values and no options to provide a catchment mask. In our case, we only want to perform simulations within a given catchment, which usually has a very irregular (not rectangular) shape. I was wondering if it would be possible to bring some (small?) modifications to TopoPyScale such that it can apply the TopoSUB clustering on an area of irregular shape.

2. From my understanding, the clustering done by TopoSUB can only deal with continuous variables (e.g. elevation, slope, etc..) as predictors. In my work, it would be greatly useful to be able to also add discrete (categorical) variables in the clustering. More specifically, the land cover type is a very important model parameter (e.g. grass, bare soil, glacier, etc..), such that I would not want grid cells with different land-cover types to belong to a same cluster. For now, I am first doing the clustering with TopoSUB, then I am adding a post-processing step (in Matlab) where I intersect the land cover map with the cluster map, which increases the number of clusters by a factor of 2 to 4. See an example in Figure 4.

3. As the work of my colleagues and myself has a clear focus on the glacier areas, we might need to increase the number of clusters for a given land-cover type (e.g. glacier). I was wondering if this is something that could also be done with TopoPyScale. See an example in Figure 5 (which was also obtained after some post-processing done in Matlab).

Sadly, I don’t have much experience with Python to directly suggest some modifications to the code to bring the changes that I would need, but I would be glad to hear your opinion on all of this. I already had some discussions with Joel Fiddes who was quite enthusiastic about this and recommended that I post this message here.

Best regards,

Achille

[ArcticSnow]

Hello @achillejouberton , welcome onboard!

Very neat work! and very nice figures!

1. yes, that would be nice to have. Currently a work around is to add an extra feature in the clustering, corresponding to a mask (0/1). This should split the clusters in two. for current work I am doing with forest, I implemented a PCA method prior to the clustering to reduce dimensionality of features, useful when using many more features against which to cluster, and in particualr if some features have colinearity. This is not yet implemented in the main TopoPyScale. So short answer: yes that would be neat, there is a work around
2. the clustering algorithm takes categorical variable in the form of translating each categories to an integer. for example: glacier=1, snow=2, rock=3, vegetation=4. This can then be thrown in the K-means algo without problem.
3. you can easily crank up the number of clusters. Actually there is a routine using the elbow method to identify a good number of clusters for your project. Look at search_number_of_clusters(). With high number of clusters, use the minibatchkmeans method that parallelize the clustering.

I hope this help you.

[ArcticSnow]

FYI, the clustering algorithm is the K-means routines from scikit-learn. Others could easily be implemented as scikit-leanr provides a common syntax for a broad range of ML algorithm.

[achillejouberton]

Thanks a lot, yes it is very helpful.

Sorry for the basic question, but how can I provide categorical variables to the clustering algorithm ?

What I have tried:

veg_file = fr'I:\Achille\T&C\Parlung4_VEG_CODE_100m_32647.tif'  # array of numbers, e.g. 3 for grass, 7 for rock, 8 for glaciers, etc..

VEG_CODE = rio.open_rasterio(veg_file).to_dataset('band')
VEG_CODE = VEG_CODE.rename({1: 'land_cover'})
VEG_data = VEG_CODE['land_cover'].to_pandas().astype('category')

ds = mp.toposub.ds_param
ds['land_cover'] = (["y", "x"], VEG_data)
ds['land_cover'] = ds['land_cover'].astype('int32')

I have tried like this and I still have clusters that are spanning several land-cover classes. I guess the data needs to be given as categorical and not 'int32' but when I tried to put " astype('category') " it led to an error.

Thanks!
Achille

[ArcticSnow]

Assign a number (int) to each category. Then you need to specify in the config file you gonna use the feature land_cover for clustering, in the dictionary. You may put feature importance of elevation to 1, rather than 4.
This should do the trick. you may also increase the feature importance of the land_cover

let me know if that works

[ArcticSnow]

oh and append the feature to mp.toposub.ds_param.

[achillejouberton]

Hello!
Here is the state of my progress:
I managed to assign a number (int) to each category, added the catchment mask and land cover features to mp.toposub.ds_param and also to the dictionary in the config file. I saw a clear difference in the clustering output, such that I am confident that the newly added features were well considered during the k-mean clustering. I also saw, though, that the centroids could get values different from the land-cover integers (for example 3.064, while land-cover categories were assigned integer values from 1 to 9). This means that some clusters are covering several land-cover types. After thinking, if I don't want a cluster to contain several land-cover classes, then it doesn't make sense to provide the land-cover feature to the k-mean algorithm, and I should rather run the k-mean algorithm directly over a given land-cover type. I made 2 tests which it turned out quite successful:

1. Using the mask feature (which is basically an array of 0 and 1 as you suggested earlier), I removed the rows of the dataframe for which the mask values are 0, in topoclass.py , at the beginning of the def extract_topo_cluster_param(self): section.

   df_param = ts.ds_to_indexed_dataframe(self.toposub.ds_param)
   index_to_drop = df_param[ df_param['mask'] == 0].index 
   index_to_keep = df_param[df_param['mask'] == 1].index
   df_param.drop(index_to_drop, inplace = True) # Drop pixels outside of mask 

The K-mean clustering is then run other this sub-sample of points, and then I adjusted the creation of the cluster map with "NaN" outside the mask:

  cluster_map = np.full(self.toposub.ds_param.slope.shape, np.nan).flatten()
  cluster_map[index_to_keep] = df_param.cluster_labels.values
  self.toposub.ds_param['cluster_labels'] = (
  ["y", "x"], np.reshape(cluster_map, self.toposub.ds_param.slope.shape))

See the figure below, it seems to work quite well, and would to solve the first issue I raised at the start of the discussion.

Then I also attempted to do that for a land-cover class, removing from the dataframe all rows for which the land-cover is not equal to - for example - 8 (with 8 = glacier).

df_param = ts.ds_to_indexed_dataframe(self.toposub.ds_param)
index_to_drop = df_param[ df_param['land_cover'] != 8].index 
index_to_keep = df_param[ df_param['land_cover'] == 8].index
df_param.drop(index_to_drop, inplace = True) # Drop pixels outside of a land-cover

It works as well, cf. figure below. Now it seems to me that to achieve what I want, I would need to loop over the land-cover class and apply the k-mean clustering for each sub-samples, and then combine the clusters in one final map. The total number of clusters (n_clusters) chosen in the config file could be distributed somehow to each land-cover types (based on their relative area?), and it could be also possible to increase the number of cluster for a given land-cover class (for e..g., glacier).
That looks pretty achievable, and my colleague Pascal Wyss who's more familiar than me with TopoPyScale and Python could give me a hand.
What are your thoughts on this Simon ?

Best,
Achille

[ArcticSnow]

you may have to 1) fork the project, 2) do changes in your fork, and 3) do a PR. Otherwise, Pascal has rights to push to the main repo.

[paswyss]

I was able to implement it :)
See commit 81298bf

(-> there are two optional configuration parameters for toposub: a mask_file and a group_file)
I will adapt the documentation later and then close this discussion ok?

[ArcticSnow]

@paswyss Very nice! Have you run the 3 examples with the new changes to check it all goes smoothly? This new feature will be super useful! Thank you

[paswyss]

Yes I have ran it.

[ArcticSnow]

then I let you close the discussion when ready. :)