fixed bug in hwsd2 convert.py and created simple visualize_netcdf.py script

Author: msteckle

HWSD2/convert.py

@@ -27,10 +27,10 @@
 #####################################################
 
 # main parameters
-# VAR = "cSoil"
-VAR = "cSoilAbove1m"
-# LAYERS = ["D1", "D2", "D3", "D4", "D5", "D6", "D7"]  # cSoil
-LAYERS = ["D1", "D2", "D3", "D4", "D5"]  # cSoilAbove1m
+VAR = "cSoil"
+# VAR = "cSoilAbove1m"
+LAYERS = ["D1", "D2", "D3", "D4", "D5", "D6", "D7"]  # cSoil
+# LAYERS = ["D1", "D2", "D3", "D4", "D5"]  # cSoilAbove1m
 POOLS = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]  # soil types
 SDATE = datetime.datetime(1960, 1, 1)
 EDATE = datetime.datetime(2022, 1, 1)
@@ -303,11 +303,13 @@ def create_netcdf(
     )
 
     # export as netcdf
-    time_range = f"{ds['time'].min().dt.year:d}{ds['time'].min().dt.month:02d}"
-    time_range += f"-{ds['time'].max().dt.year:d}{ds['time'].max().dt.month:02d}"
     ds.to_netcdf(
-        "{variable}_{frequency}_{source_id}_{time_mark}.nc".format(
-            variable=var, frequency="fx", source_id="HWSD2", time_mark=time_range
+        "{variable}_{frequency}_{source_id}_{st_date}-{en_date}.nc".format(
+            variable=var,
+            frequency="fx",
+            source_id="HWSD2",
+            st_date=sdate.strftime("%Y%m%d"),
+            en_date=edate.strftime("%Y%m%d"),
         ),
         encoding={VAR: {"zlib": True}},
     )

scripts/visualize_netcdf.py

@@ -0,0 +1,48 @@
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import xarray as xr
+
+# set parameters
+wdir = "../HWSD2"
+file_path = f"{wdir}/cSoil_fx_HWSD2_19600101-20220101.nc"
+tstep = 0
+vmin = 0  # Minimum value for colormap
+vmax = 100  # Maximum value for colormap
+var = "cSoil"
+
+# open netcdf and select time step
+base_name = Path(file_path).stem
+data = xr.open_dataset(file_path)
+da = data[var].isel(time=tstep)
+
+# create and save full map
+plt.figure(figsize=(10, 6))
+p = da.plot(vmin=vmin, vmax=vmax)
+plt.savefig(f"{wdir}/{base_name}_timestep_{tstep}.png", dpi=300, bbox_inches="tight")
+plt.close()
+
+#### Create zoomed-in map ####
+# Define Southeastern US bounding box
+lon_min, lon_max = -95, -75
+lat_min, lat_max = 25, 37
+
+# Determine proper slicing directions
+lat_vals = da["lat"].values
+lon_vals = da["lon"].values
+
+lat_slice = (
+    slice(lat_min, lat_max) if lat_vals[0] < lat_vals[-1] else slice(lat_max, lat_min)
+)
+lon_slice = (
+    slice(lon_min, lon_max) if lon_vals[0] < lon_vals[-1] else slice(lon_max, lon_min)
+)
+
+# Clip the data dynamically
+da_se = da.sel(lon=lon_slice, lat=lat_slice)
+
+# Plot the clipped region
+plt.figure(figsize=(8, 6))
+p = da_se.plot(vmin=vmin, vmax=vmax)
+plt.savefig(f"{wdir}/{base_name}_SE-US_zoom.png", dpi=300, bbox_inches="tight")
+plt.close()