Habitat suitability for the Great Basin Bristlecone Pine#331
Conversation
|
@MarJess @rainihelmstadter you are assigned to review this |
| # - Soil pH: >7.0 | ||
| # - Elevation: 8,000 - 12,000 ft | ||
| # - Aspect: South or west preferred (values range from 0.7 to 1.0) | ||
| # - Slope: 10 - 50 degrees | ||
| # - Precipitation: 30-70 cm per year | ||
| # - Max average annual temperature: 10.0 through 18.0 degrees C | ||
| # - Min average annual temperature -4.0 through 4.0 degrees C |
There was a problem hiding this comment.
This is a very nice summary for later when it comes to the fuzzy logic model.
Very clear and clean.
| ### read in the file | ||
| pa_shp = gpd.read_file(pa_path, layer = "PADUS4_1Fee_State_NV") | ||
|
|
||
| ### check the crs | ||
| pa_shp.crs | ||
|
|
||
| ### assign crs to match GBIF | ||
| pa_shp = pa_shp.to_crs(epsg = 4326) | ||
|
|
||
| ### Fix invalid geometries to prevent spatial operation errors | ||
| pa_shp['geometry'] = pa_shp['geometry'].apply( | ||
| lambda geom: geom.make_valid() if not isinstance(geom, | ||
| MultiPolygon) and not geom.is_valid else geom) | ||
|
|
||
| ### drop remaining invalid geoms | ||
| pa_shp = pa_shp[pa_shp.geometry.is_valid] | ||
|
|
||
| ### drop missing geometry | ||
| pa_shp = pa_shp.dropna(subset=['geometry']) | ||
|
|
||
| ### subset the data | ||
| subset_pa = pa_shp | ||
|
|
There was a problem hiding this comment.
You could wrap these steps in a function, though in our case this is only done once, so it might not be necessary.
| def download_polaris(site_name, site_gdf, soil_prop, stat, soil_depth, | ||
| plot_path, plot_title, data_dir, plots_dir): |
There was a problem hiding this comment.
Thats a nice solution to have a wrapper function which executes the other functions. This makes it clear and not like a spaghetti code.
Optional, this could also be replaced by a POLARIS class with the functions as methods.
|
|
||
| # %% | ||
| ### function to get aspect layer from elevation data, and plot it | ||
| def plot_aspect(site_name, srtm_da, site_gdf, topo_site_dir): |
There was a problem hiding this comment.
a little bit misleading function name. At first I would assume that this function is only for plotting.
But I actually gets the aspect and plots it.
| # %% | ||
| ### function to get the slope | ||
|
|
||
| def plot_slope(site_name, srtm_da, site_gdf, topo_site_dir): |
| ) | ||
|
|
||
|
|
||
| def check_alignment(rasters): |
There was a problem hiding this comment.
thats a neat little helper function.
| # - Cold and Dry: GFDL-ESM2M | ||
| # - Cold and Wet: CNRM-CM5 | ||
|
|
||
| # %% [markdown] |
There was a problem hiding this comment.
The following is a very nice and detail results part.
Thanks for that!
MarJess
left a comment
There was a problem hiding this comment.
Thank you max for this very nice analysis.
It is well structured. There are a few steps, which may can be included in a loop to make the code more compact. But honestly, this sometimes makes it harder to read.
ryme1295
left a comment
There was a problem hiding this comment.
I think this was great and an example I'll try to follow during my final project!
| # | ||
| # ### Study Sites | ||
| # For this study, I will select two sites where Great Basin Bristlecone Pines are common. The first site is [Great Basin National Park](https://www.nps.gov/grba/index.htm) in the Northeastern part of Nevada. The other site is a portion of [Humboldt-Toiyabe National Forest](https://www.nationalforests.org/forest/humboldt-toiyabe-national-forest/) just west of Las Vegas. This will allow for comparison of how habitat suitability is changing in the northern and southern parts of the state for this species. | ||
| # |
There was a problem hiding this comment.
Something I am working on doing better is a summary of what I am using and the goal in mind. This could also turn into an area where you mention what could give you difficulty.
| fill_color = None, | ||
| line_color = "white", | ||
| frame_width = 600 | ||
| ) |
There was a problem hiding this comment.
Here, I liked the step of using satellite imagery, adding a title, and drawing shapes with white outlines. You could add a legend or labels so users understand what the outlines represent.
| aspect = info["aspect"] | ||
| soil = info["soil"] | ||
| boundary = info["boundary"] | ||
|
|
There was a problem hiding this comment.
Nice job gathering the variables; it makes your process easy to see.
optional- include precipitation.
|
|
||
| # %% | ||
| ### plot it | ||
| gbif_gdf.hvplot( |
There was a problem hiding this comment.
This plot shows the entire globe, vs just the areas I think you're interested in, likely due to missing data in the occurrence dataset. The telltale sign is that there's an occurrence near the west coast of Africa, which is where (0,0) plots. I'd suggest filtering the gbif_gdf before plotting to remove occurrences without location data.
| ) | ||
|
|
||
| # %% | ||
| ### remove small polygons inside boundary to prevent errors |
There was a problem hiding this comment.
Nice work isolating just the area of HTNF you were interested in, and good catch on removing the small polygons!
| buffer (float): buffer around bounding box (degrees) | ||
|
|
||
| Returns: | ||
| tuple: buffered bounding box |
There was a problem hiding this comment.
Nice thinking to export the buffered bounding box so you can use it later. Super small note, but I think these returns could be named something more applicable, like bounding_box (tuple) and file_pattern (str) or something.
| fig, ax = plt.subplots(figsize=(8,6)) | ||
|
|
||
| ### plot raster | ||
| aspect_plot = aspect_da.plot(ax=ax, cmap='terrain') |
There was a problem hiding this comment.
You could consider using a colormap that is cyclic here, as aspects at the ends of the range are actually convergent. Some nice documentation here (if that's helpful): https://matplotlib.org/stable/users/explain/colors/colormaps.html
| step (int): number of years per range (default is 5) | ||
|
|
||
| Returns: | ||
| list of str: list of year range strings in "startYear_endYear" format |
There was a problem hiding this comment.
Just a thought here again, you could state the name of the variable you're returning, as well as the datatype. I might be stylistically incorrect here but I feel like it would make the docstring a hair more informative.
rainihelmstadter
left a comment
There was a problem hiding this comment.
Nice work Max! Your code is overall clear, organized, and well-documented. I liked how you added the check functions and brought the GBIF occurrence data back at the end. I'd recommend just a bit of editing to ensure clarity throughout - see my comments on the GBIF occurrence plot and returns sections of some of your docstrings.
| print(f"\nAlignment check for {site_name}") | ||
| first_model = models[0] | ||
|
|
||
| check_alignment([ |
There was a problem hiding this comment.
Nice check function on alignment!
| Returns: | ||
| None | ||
|
|
||
| Notes: |
There was a problem hiding this comment.
I appreciate the notes you've left on a number of these functions! You preserve readability while adding context in a nice way.
| plt.show() | ||
|
|
||
| # %% [markdown] | ||
| # ### Plot the final habitat suitability maps with GBIF occurrence data |
There was a problem hiding this comment.
I liked that you included the GBIF occurrence data! However, the tight clustering of some of the occurrences, especially in HTNF, made it hard to see how suitability will change in that area because the occurrences took up a lot of visual space. Maybe there's a way to show the occurrences at the beginning of the plots, or change the color/alpha to ensure that the data underneath isn't changed?
Here is a link to the html version to view plots and outputs without having to run the code yourself.
https://maxwarnock.github.io/projects/habitat_suitability.html