@@ -8,6 +8,7 @@ description: >
88 In this tutorial we will run overland flow
99 and erosion simulation (SIMWE) split by subwatersheds
1010 in parallel.
11+ thumbnail : SIMWE_images/thumbnail.webp
1112format :
1213 html :
1314 toc : true
@@ -677,8 +678,7 @@ import grass.script as gs
677678
678679
679680def compute(huc12):
680- simulation_time = 300
681- print(huc12)
681+ simulation_time = 180
682682 gs.run_command(
683683 "v.extract",
684684 input="river_basins",
@@ -698,6 +698,7 @@ def compute(huc12):
698698 output_directory=os.getcwd(),
699699 title_filter="20240510",
700700 output_name=f"ned_{huc12}",
701+ quiet=True,
701702 env=env,
702703 )
703704 with gs.RegionManager(raster=f"ned_{huc12}", env=env):
@@ -778,15 +779,16 @@ def compute(huc12):
778779 )
779780 # Erosion: extract negative values and convert to positive mass [kg/m^2]
780781 gs.mapcalc(
781- f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {simulation_time * 60}, null() )",
782+ f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {simulation_time * 60}, 0 )",
782783 env=env,
783784 )
784785 erosion = gs.parse_command(
785- "r.univar", map=f"erosion_{huc12}", format="json", env=env
786+ "r.univar", map=f"erosion_{huc12}", format="json", flags="e", env=env
786787 )
787788 return {
788789 "huc12": huc12,
789790 "erosion_mean": erosion["mean"],
791+ "erosion_median": erosion["median"],
790792 "erosion_total": erosion["sum"] * region["nsres"] * region["ewres"],
791793 }
792794
@@ -796,13 +798,11 @@ if __name__ == "__main__":
796798 "records"
797799 ]
798800 huc12s = [basin["huc12"] for basin in basins]
799-
800801 # set the number of processes:
801802 with Pool(processes=cpu_count()) as pool:
802803 result = list(tqdm(pool.imap(compute, huc12s), total=len(huc12s)))
803804 with open("result.json", "w") as fp:
804805 json.dump(result, fp)
805-
806806```
807807
808808Now execute the script, it will take some time.
@@ -815,18 +815,14 @@ Now execute the script, it will take some time.
815815Load the resulting JSON file into a pandas dataframe and normalize the results, so that we can easily compare the subwatersheds.
816816
817817``` {python}
818+ import pandas as pd
818819import json
819820
820821with open("result.json") as f:
821822 stats = json.load(f)
822823
823824df = pd.DataFrame(stats)
824- df["normalized_erosion"] = df.erosion_mean / max(
825- max(df.erosion_mean), max(df.deposition_mean)
826- )
827- df["normalized_deposition"] = df.deposition_mean / max(
828- max(df.erosion_mean), max(df.deposition_mean)
829- )
825+ df["normalized_erosion"] = df.erosion_mean / max(df.erosion_mean)
830826df
831827```
832828
@@ -851,15 +847,16 @@ Let's visualize the data with folium:
851847import folium
852848import branca.colormap as cm
853849
854- metric = "normalized_deposition "
850+ metric = "normalized_erosion "
855851color_scale = cm.linear.OrRd_09.scale(min(gdf[metric]), max(gdf[metric]))
856- color_scale.caption = "Min-max normalized erosion"
852+ color_scale.caption = "Min-max normalized mean erosion"
853+ # get region center as latitude, longitude
857854region = gs.parse_command("g.region", vector="river_basins", format="json", flags="b")
858855m = folium.Map(location=[region["ll_clat"], region["ll_clon"]], zoom_start=9)
859856
860857tooltip = folium.GeoJsonTooltip(
861- fields=["huc12", "normalized_erosion", "normalized_deposition" ],
862- aliases=["HUC12", "Normalized erosion:", "Normalized deposition :"],
858+ fields=["huc12", metric ],
859+ aliases=["HUC12", "Normalized mean erosion :"],
863860 localize=True,
864861 sticky=False,
865862 labels=True,
@@ -876,9 +873,8 @@ g = folium.GeoJson(
876873 gdf,
877874 style_function=lambda x: {
878875 "fillColor": color_scale(x["properties"][metric]),
879- "color": "black",
880- "weight": 1,
881- "fillOpacity": 0.4,
876+ "color": "black", "weight": 1,
877+ "fillOpacity": 0.5,
882878 },
883879 tooltip=tooltip,
884880).add_to(m)
@@ -887,5 +883,5 @@ folium.LayerControl().add_to(m)
887883m
888884```
889885
890-
886+ ![ ] ( SIMWE_images/webmap.webp )
891887
0 commit comments