From 44eb56c03d9a89315e3c0e998bc3dd8e007282e8 Mon Sep 17 00:00:00 2001 From: Anthony Date: Tue, 27 Aug 2024 08:42:59 -0400 Subject: [PATCH 1/2] Added land analysis test scripts --- .../freanalysis_aerosol/__init__.py | 3 +- freanalysis_land/README.md | 1 + freanalysis_land/freanalysis_land/__init__.py | 0 freanalysis_land/freanalysis_land/land.py | 243 ++++++++++++++++++ freanalysis_land/pyproject.toml | 35 +++ tests/my_land_analysis.py | 5 + 6 files changed, 286 insertions(+), 1 deletion(-) create mode 100644 freanalysis_land/README.md create mode 100644 freanalysis_land/freanalysis_land/__init__.py create mode 100644 freanalysis_land/freanalysis_land/land.py create mode 100644 freanalysis_land/pyproject.toml create mode 100644 tests/my_land_analysis.py diff --git a/freanalysis_aerosol/freanalysis_aerosol/__init__.py b/freanalysis_aerosol/freanalysis_aerosol/__init__.py index bc8971a..734a3e6 100644 --- a/freanalysis_aerosol/freanalysis_aerosol/__init__.py +++ b/freanalysis_aerosol/freanalysis_aerosol/__init__.py @@ -11,7 +11,7 @@ @dataclass class Metadata: """Helper class that stores the metadata needed by the plugin.""" - frequency: str = "monthly" + frequency: str = "monthly"object realm: str = "atmos" @staticmethod @@ -207,3 +207,4 @@ def run_analysis(self, catalog, png_dir, reference_catalog=None, config={}): figure.save(Path(png_dir) / f"{name}.png") figure_paths.append(Path(png_dir)/ f"{name}.png") return figure_paths +object \ No newline at end of file diff --git a/freanalysis_land/README.md b/freanalysis_land/README.md new file mode 100644 index 0000000..ca4e23a --- /dev/null +++ b/freanalysis_land/README.md @@ -0,0 +1 @@ +This is the README. \ No newline at end of file diff --git a/freanalysis_land/freanalysis_land/__init__.py b/freanalysis_land/freanalysis_land/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/freanalysis_land/freanalysis_land/land.py b/freanalysis_land/freanalysis_land/land.py new file mode 100644 index 0000000..3d1f62b --- /dev/null +++ b/freanalysis_land/freanalysis_land/land.py @@ -0,0 +1,243 @@ +import json +from analysis_scripts import AnalysisScript +import intake +import matplotlib.pyplot as plt +import cartopy +import cartopy.crs as ccrs +import datetime +import pandas as pd +import re +import xarray as xr + + +class LandAnalysisScript(AnalysisScript): + """Abstract base class for analysis scripts. User-defined analysis scripts + should inhert from this class and override the requires and run_analysis methods. + + Attributes: + description: Longer form description for the analysis. + title: Title that describes the analysis. + """ + def __init__(self): + """Instantiates an object. The user should provide a description and title.""" + self.description = "This is for analysis of land model (stand-alone)" + self.title = "Soil Carbon" + + def requires(self): + """Provides metadata describing what is needed for this analysis to run. + + Returns: + A json string describing the metadata. + """ + raise NotImplementedError("you must override this function.") + return json.dumps("{json of metadata MDTF format.}") + + def global_map(self,dataset,var,dates,plt_time=None,colormap='viridis',title=''): + """ + Generate a global map and regional subplots for a specified variable from an xarray dataset. + + This function creates a global map and several regional subplots (North America, South America, Europe, Africa, Asia, and Australia) + using the specified variable from an xarray dataset. The generated map will be saved as a PNG file. + + Parameters: + ---------- + dataset : xarray.Dataset + The input xarray dataset containing the variable to be plotted. + + var : str + The name of the variable in the dataset to be plotted. + + dates: list + The list of dates from the dataframe, converted to period index. + + plt_time : int, optional + The time index to plot from the variable data. Defaults to the length of `dates` - 1, or last date in dataset. + + colormap : str, optional + The colormap to use for plotting the data. Defaults to 'viridis'. + + title : str, optional + The title for the figure. Defaults to an empty string. + + Returns: + ------- + fig : matplotlib figure object + The function returns the plot for saving + + Notes: + ----- + The function uses Cartopy for map projections and Matplotlib for plotting. + Ensure Cartopy and Matplotlib are installed in your Python environment. + + The output file is saved with the format '_global_map.png'. + + Examples: + -------- + global_map(my_dataset, 'mrso', plt_time=0, colormap='coolwarm', title='Global Temperature Map') + """ + lon = dataset.lon + lat = dataset.lat + + if plt_time == None: plt_time=len(dates)-1 + + data = dataset[var][plt_time].values + projection = ccrs.PlateCarree() + fig = plt.figure(figsize=(8.5, 11)) + + # Global map + ax_global = fig.add_subplot(3, 1, 1, projection=projection) + ax_global.set_title('Global Map') + mesh = ax_global.pcolormesh(lon, lat, data, transform=projection, cmap=colormap) + ax_global.coastlines() + + # List of bounding boxes for different continents (min_lon, max_lon, min_lat, max_lat) + regions = { + 'North America': [-170, -47, 0, 85], + 'South America': [-90, -30, -60, 15], + 'Europe': [-10, 60, 30, 75], + 'Africa': [-20, 50, -35, 37], + 'Asia': [60, 150, 5, 75], + 'Australia': [110, 180, -50, 0] + } + + # Create subplots for each region + for i, (region, bbox) in enumerate(regions.items(), start=1): + ax = fig.add_subplot(3, 3, i + 3, projection=projection) + ax.set_extent(bbox, crs=projection) + ax.set_title(region) + ax.pcolormesh(lon, lat, data, transform=projection, cmap=colormap) + ax.coastlines() + ax.add_feature(cartopy.feature.BORDERS) + + # Add colorbar + fig.colorbar(mesh, ax=ax_global, orientation='horizontal', pad=0.05, aspect=50) + plt.suptitle(title) + plt.tight_layout() + + return fig + # plt.savefig(var+'_global_map.png') + # plt.close() + + def timeseries(self, dataset,var,dates_period,var_range=None,minlon = 0,maxlon = 360,minlat = -90,maxlat=90,timerange=None,title=''): + ''' + Generate a time series plot of the specified variable from a dataset within a given geographic and temporal range. + + + Parameters: + ----------- + dataset : xarray.Dataset + The dataset containing the variable to be plotted. + var : str + The name of the variable to plot from the dataset. + dates_period : pandas.DatetimeIndex + The dates for the time series data. + var_range : tuple of float, optional + The range of variable values to include in the plot (min, max). If not provided, the default range is (0, inf). + minlon : float, optional + The minimum longitude to include in the plot. Default is 0. + maxlon : float, optional + The maximum longitude to include in the plot. Default is 360. + minlat : float, optional + The minimum latitude to include in the plot. Default is -90. + maxlat : float, optional + The maximum latitude to include in the plot. Default is 90. + timerange : tuple of int, optional + The range of years to plot (start_year, end_year). If not provided, all available years in the dataset will be plotted. + title : str, optional + The title of the plot. Default is an empty string. + + Returns: + -------- + matplotlib.figure.Figure + The figure object containing the generated time series plot. + + Notes: + ------ + The function filters the dataset based on the provided variable range, longitude, and latitude bounds. It then + calculates the monthly and annual means of the specified variable and plots the seasonal and annual means. + + ''' + if var_range != None: + data_filtered = dataset.where((dataset[var] > var_range[0]) & (dataset[var] <= var_range[1]) & + (dataset.lat >= minlat) & (dataset.lon >= minlon) & + (dataset.lat <= maxlat) & (dataset.lon <= maxlon)) + else: + data_filtered = dataset.where((dataset[var] > 0) & + (dataset.lat >= minlat) & (dataset.lon >= minlon) & + (dataset.lat <= maxlat) & (dataset.lon <= maxlon)) + data_filtered['time'] = dates_period + + data_df = pd.DataFrame(index = dates_period) + data_df['monthly_mean'] = data_filtered.resample(time='YE').mean(dim=['lat','lon'],skipna=True)[var].values + data_df['monthly_shift'] = data_df['monthly_mean'].shift(1) + + if timerange != None: + ys, ye = (str(timerange[0]),str(timerange[1])) + plot_df = data_df.loc[f'{ys}-1-1':f'{ye}-1-1'] + else: + plot_df = data_df + + fig, ax = plt.subplots() + plot_df.resample('Q').mean()['monthly_shift'].plot(ax=ax,label='Seasonal Mean') + plot_df.resample('Y').mean()['monthly_mean'].plot(ax=ax,label='Annual Mean') + plt.legend() + plt.title(title) + plt.xlabel('Years') + return fig + + def run_analysis(self, catalog, png_dir, reference_catalog=None): + """Runs the analysis and generates all plots and associated datasets. + + Args: + catalog: Path to a model output catalog. + png_dir: Directory to store output png figures in. + reference_catalog: Path to a catalog of reference data. + + Returns: + A list of png figures. + """ + print ('WARNING: THESE FIGURES ARE FOR TESTING THE NEW ANALYSIS WORKFLOW ONLY AND SHOULD NOT BE USED IN ANY OFFICIAL MANNER FOR ANALYSIS OF LAND MODEL OUTPUT.') + col = intake.open_esm_datastore(catalog) + df = col.df + + # Soil Carbon + var = 'cSoil' + print ('Soil Carbon Analysis') + cat = col.search(variable_id=var,realm='land_cmip') + other_dict = cat.to_dataset_dict(cdf_kwargs={'chunks':{'time':12},'decode_times':False}) + combined_dataset = xr.concat(list(dict(sorted(other_dict.items())).values()), dim='time') + + # Other data: + # land_static_file = re.search('land_static:\s([\w.]*)',combined_dataset.associated_files).group(1) + # STATIC FILES SHOULD BE PART OF THE CATALOG FOR EASY ACCESS + + # Select Data and plot + dates = [datetime.date(1,1,1) + datetime.timedelta(d) for d in combined_dataset['time'].values] # Needs to be made dynamic + dates_period = pd.PeriodIndex(dates,freq='D') + + sm_fig = self.global_map(combined_dataset,var,dates,title='Soil Carbon Content (kg/m^2)') + plt.savefig(png_dir+var+'_global_map.png') + plt.close() + + ts_fig = self.timeseries(combined_dataset,var,dates_period,title='Global Average Soil Carbon') + plt.savefig(png_dir+var+'_global_ts.png') + plt.close() + + # Soil Moisture + var = 'mrso' + print ('Soil Moisture Analysis') + cat = col.search(variable_id=var,realm='land_cmip') + other_dict = cat.to_dataset_dict(cdf_kwargs={'chunks':{'time':12},'decode_times':False}) + combined_dataset = xr.concat(list(dict(sorted(other_dict.items())).values()), dim='time') + + # Other data: + # soil_area_file = re.search('soil_area:\s([\w.]*)',combined_dataset.associated_files).group(1) + # STATIC FILES SHOULD BE PART OF THE CATALOG FOR EASY ACCESS + + # Select Data and plot + dates = [datetime.date(1,1,1) + datetime.timedelta(d) for d in combined_dataset['time'].values] # Needs to be made dynamic + dates_period = pd.PeriodIndex(dates,freq='D') + + sm_fig = self.global_map(combined_dataset,var,dates,title='Soil Moisture (kg/m^2)') + plt.savefig(png_dir+var+'_global_map.png') + plt.close() diff --git a/freanalysis_land/pyproject.toml b/freanalysis_land/pyproject.toml new file mode 100644 index 0000000..c0cc45d --- /dev/null +++ b/freanalysis_land/pyproject.toml @@ -0,0 +1,35 @@ +[build-system] +requires = [ + "setuptools >= 40.9.0", +] +build-backend = "setuptools.build_meta" + +[project] +name = "freanalysis_land" +version = "0.1" +dependencies = [ + "setuptools", + "intake", + "intake-esm", + "xarray", + "matplotlib", + "cartopy", + "pandas", + "xarray" + +] +requires-python = ">= 3.6" +authors = [ + {name = "developers"}, +] +maintainers = [ + {name = "Anthony Preucil", email = "Anthony.Preucil@noaa.gov"}, +] +description = "Land Analysis for GFDL stand-alone model" +readme = "README.md" +classifiers = [ + "Programming Language :: Python" +] + +[project.urls] +repository = "https://github.com/NOAA-GFDL/analysis-scripts.git" diff --git a/tests/my_land_analysis.py b/tests/my_land_analysis.py new file mode 100644 index 0000000..ed4104b --- /dev/null +++ b/tests/my_land_analysis.py @@ -0,0 +1,5 @@ + + +from freanalysis_land.land import LandAnalysisScript +land = LandAnalysisScript() +land.run_analysis("/work/a2p/lm4p2sc_catalog.json","/work/a2p/") \ No newline at end of file From 7f018196fb79d08efb364e03df7f571884c01dff Mon Sep 17 00:00:00 2001 From: Anthony Date: Mon, 16 Sep 2024 15:25:53 -0400 Subject: [PATCH 2/2] Responding to suggestion from Ian on my pull request. --- freanalysis_aerosol/freanalysis_aerosol/__init__.py | 5 ++--- freanalysis_land/freanalysis_land/land.py | 10 +++++----- tests/my_land_analysis.py | 5 ----- tests/test_freanalysis_land.py | 10 ++++++++++ 4 files changed, 17 insertions(+), 13 deletions(-) delete mode 100644 tests/my_land_analysis.py create mode 100644 tests/test_freanalysis_land.py diff --git a/freanalysis_aerosol/freanalysis_aerosol/__init__.py b/freanalysis_aerosol/freanalysis_aerosol/__init__.py index 734a3e6..613ed76 100644 --- a/freanalysis_aerosol/freanalysis_aerosol/__init__.py +++ b/freanalysis_aerosol/freanalysis_aerosol/__init__.py @@ -11,7 +11,7 @@ @dataclass class Metadata: """Helper class that stores the metadata needed by the plugin.""" - frequency: str = "monthly"object + frequency: str = "monthly" realm: str = "atmos" @staticmethod @@ -206,5 +206,4 @@ def run_analysis(self, catalog, png_dir, reference_catalog=None, config={}): ) figure.save(Path(png_dir) / f"{name}.png") figure_paths.append(Path(png_dir)/ f"{name}.png") - return figure_paths -object \ No newline at end of file + return figure_paths \ No newline at end of file diff --git a/freanalysis_land/freanalysis_land/land.py b/freanalysis_land/freanalysis_land/land.py index 3d1f62b..0bd3f0a 100644 --- a/freanalysis_land/freanalysis_land/land.py +++ b/freanalysis_land/freanalysis_land/land.py @@ -11,8 +11,8 @@ class LandAnalysisScript(AnalysisScript): - """Abstract base class for analysis scripts. User-defined analysis scripts - should inhert from this class and override the requires and run_analysis methods. + """A class for performing various analysis tasks relating to GFDL land model output, + inherits from the AnalysisScipt base class. Attributes: description: Longer form description for the analysis. @@ -78,7 +78,7 @@ def global_map(self,dataset,var,dates,plt_time=None,colormap='viridis',title='') lon = dataset.lon lat = dataset.lat - if plt_time == None: plt_time=len(dates)-1 + if plt_time is None: plt_time=len(dates)-1 data = dataset[var][plt_time].values projection = ccrs.PlateCarree() @@ -157,7 +157,7 @@ def timeseries(self, dataset,var,dates_period,var_range=None,minlon = 0,maxlon = calculates the monthly and annual means of the specified variable and plots the seasonal and annual means. ''' - if var_range != None: + if var_range is not None: data_filtered = dataset.where((dataset[var] > var_range[0]) & (dataset[var] <= var_range[1]) & (dataset.lat >= minlat) & (dataset.lon >= minlon) & (dataset.lat <= maxlat) & (dataset.lon <= maxlon)) @@ -171,7 +171,7 @@ def timeseries(self, dataset,var,dates_period,var_range=None,minlon = 0,maxlon = data_df['monthly_mean'] = data_filtered.resample(time='YE').mean(dim=['lat','lon'],skipna=True)[var].values data_df['monthly_shift'] = data_df['monthly_mean'].shift(1) - if timerange != None: + if timerange is not None: ys, ye = (str(timerange[0]),str(timerange[1])) plot_df = data_df.loc[f'{ys}-1-1':f'{ye}-1-1'] else: diff --git a/tests/my_land_analysis.py b/tests/my_land_analysis.py deleted file mode 100644 index ed4104b..0000000 --- a/tests/my_land_analysis.py +++ /dev/null @@ -1,5 +0,0 @@ - - -from freanalysis_land.land import LandAnalysisScript -land = LandAnalysisScript() -land.run_analysis("/work/a2p/lm4p2sc_catalog.json","/work/a2p/") \ No newline at end of file diff --git a/tests/test_freanalysis_land.py b/tests/test_freanalysis_land.py new file mode 100644 index 0000000..dfcd8ac --- /dev/null +++ b/tests/test_freanalysis_land.py @@ -0,0 +1,10 @@ + + +from freanalysis_land.land import LandAnalysisScript + +def test_land_analysis_script(): + land = LandAnalysisScript() + land.run_analysis("/work/a2p/lm4p2sc_GSWP3_hist_irr_catalog.json","/work/a2p/") + + +test_land_analysis_script() \ No newline at end of file