Plotting extent

[lgrimley]

I have a model that covers more than one UTM zone. When I use the plot_basemap() function it cuts off and only shows the area of the model that is in the UTM zone that the model as been assigned to. Is there a way for it to clip to the extent of the region instead?

[DirkEilander]

Thanks for the suggestion @lgrimley! I'll need to figure out why this happens, but agree this should be changed. If it is urgent, as a work around, you can copy the plot_basemaps method and adapt/delete the line where the plot extent is set https://github.com/Deltares/hydromt_sfincs/blob/main/hydromt_sfincs/plots.py#L189

[DirkEilander]

Currently the approach is to use the projection of the data. This could lead to negative coordinates if a model spans multiple UTM zones with cartopy doesn't seem to like (this assumption needs to be checked). Instead of a UTM zone it would perhaps be better if the model uses another CRS that can describe the entire model in that case.

@lgrimley can you check if this error persists if you reproject the model grid to e.g. WGS84 and then plot it?

sf = SfincsModel(<root>, mode='r')
sf._grid = sf.grid.raster.reproject('epsg:4326')
sf.plot_basemap()

[roeldegoede]

Finally got to work on this issue again, and found some things worth mentioning:

• The solution proposed by Dirk works, but only if you also reproject the geoms and forcing (used for boundary points) OR use plot_geoms=False. Otherwise, those will remain in the old CRS and plotting results in vague errors of a very large figure.
• We used the cartopy ccrs.epsg() methods, which have a very narrow range of application for UTM zones, whereas the ccrs.UTM() method does allow for a much larger margin around the UTM zone. This I will implement in the plotting script.
• I'll also implement a warning for when your model domain falls outside of the extent of the CRS.

Something to consider adding are the following lines of code in plots.py, which would automaticcaly reproject the model to WGS 84 in case the model exceeds the CRS extent:

# reproject ds to EPSG:4326
ds = ds.raster.reproject('epsg:4326')
# reproject all geoms to EPSG:4326
for name, gdf in geoms.items():
   geoms[name] = gdf.to_crs("epsg:4326")
# update proj_crs, proj_str
proj_crs = CRS.from_epsg(4326)
proj_str = proj_crs.name
#(also update bounds and extent!)
bounds = 
extent = 

[lgrimley]

Thanks everyone! Somehow I missed all of Dirk's messaged and am just now seeing them (oops). I can give these suggestions a try

[roeldegoede]

Perfect! I think with the merge #226, plotting of these models should have improved already. If it still doesn't fit in the figure, it should provide a more clear warning now that it doesn't fit in this UTM zone, and that you should reproject it to WGS84 for plotting purposes (with the code provided above)