As far as I can see in the codebase, you're not taking care of this issue neither in the reproject_raster nor in the download_elevation, though.
What I meant was that the output of download_elevation is a dem at 4326 (non-project CRS) while all your other computations are in UTM (a projected CRS). So, at the end of the day, you are reprojecting the downloaded dem to UTM which ends up producing artificial nulls. My workaround is to get the dem data for a buffered bbox so once you reproject and clip it to the original bounds, the artificial nulls get clipped and the final reprojected dem ends up not having nulls. I recommend looking at your intermedia results in the workflow to see if that's not the case.
Something like this should take care of the issue:
bbox = gdf.buffer(20*30).to_crs(4326).total_bounds
dem = download_elevation(bbox)
dem = dem.rio.reproject(gdf.crs).fillna(dem.rio.nodata)
dem = dem.rio.clip_box(*gdf.total_bounds).astype("int16")
dem.rio.to_raster("dem.tif")
Originally posted by @cheginit in https://github.com/ImperialCollegeLondon/SWMManywhere/issues/263#issuecomment-2314188321