climate-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From jo...@apache.org
Subject [1/6] git commit: Initial script creation
Date Mon, 27 Oct 2014 19:16:28 GMT
Repository: climate
Updated Branches:
  refs/heads/master 6d47a5783 -> a69f5545f


Initial script creation


Project: http://git-wip-us.apache.org/repos/asf/climate/repo
Commit: http://git-wip-us.apache.org/repos/asf/climate/commit/a0e5d0bd
Tree: http://git-wip-us.apache.org/repos/asf/climate/tree/a0e5d0bd
Diff: http://git-wip-us.apache.org/repos/asf/climate/diff/a0e5d0bd

Branch: refs/heads/master
Commit: a0e5d0bd7cfeb05fa82b725fcd5232ac2b043971
Parents: 0951a93
Author: cgoodale <sigep311@gmail.com>
Authored: Mon Jun 23 12:19:47 2014 -0700
Committer: cgoodale <sigep311@gmail.com>
Committed: Mon Jun 23 12:19:47 2014 -0700

----------------------------------------------------------------------
 examples/model_ensemble_to_rcmed.py | 185 +++++++++++++++++++++++++++++++
 1 file changed, 185 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/climate/blob/a0e5d0bd/examples/model_ensemble_to_rcmed.py
----------------------------------------------------------------------
diff --git a/examples/model_ensemble_to_rcmed.py b/examples/model_ensemble_to_rcmed.py
new file mode 100644
index 0000000..33f1219
--- /dev/null
+++ b/examples/model_ensemble_to_rcmed.py
@@ -0,0 +1,185 @@
+# Licensed to the Apache Software Foundation (ASF) under one
+# or more contributor license agreements.  See the NOTICE file
+# distributed with this work for additional information
+# regarding copyright ownership.  The ASF licenses this file
+# to you under the Apache License, Version 2.0 (the
+# "License"); you may not use this file except in compliance
+# with the License.  You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing,
+# software distributed under the License is distributed on an
+# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+# KIND, either express or implied.  See the License for the
+# specific language governing permissions and limitations
+# under the License.
+
+import datetime
+import urllib
+from os import path
+
+import numpy as np
+
+import ocw.data_source.local as local
+import ocw.data_source.rcmed as rcmed
+from ocw.dataset import Bounds as Bounds
+import ocw.dataset_processor as dsp
+import ocw.evaluation as evaluation
+import ocw.metrics as metrics
+import ocw.plotter as plotter
+
+# File URL leader
+FILE_LEADER = "http://zipper.jpl.nasa.gov/dist/"
+# This way we can easily adjust the time span of the retrievals
+YEARS = 3
+# Two Local Model Files 
+FILE_1 = "AFRICA_KNMI-RACMO2.2b_CTL_ERAINT_MM_50km_1989-2008_tasmax.nc"
+FILE_2 = "AFRICA_UC-WRF311_CTL_ERAINT_MM_50km-rg_1989-2008_tasmax.nc"
+# Filename for the output image/plot (without file extension)
+OUTPUT_PLOT = "model_ensemble_tasmax_africa_bias_monthly"
+
+# Download necessary NetCDF file if not present
+if path.exists(FILE_1):
+    pass
+else:
+    urllib.urlretrieve(FILE_LEADER + FILE_1, FILE_1)
+
+if path.exists(FILE_2):
+    pass
+else:
+    urllib.urlretrieve(FILE_LEADER + FILE_2, FILE_2)
+
+
+""" Step 1: Load Local NetCDF File into OCW Dataset Objects """
+# Load local knmi model data
+knmi_dataset = local.load_file(FILE_1, "tasmax")
+knmi_dataset.name = "AFRICA_KNMI-RACMO2.2b_CTL_ERAINT_MM_50km_1989-2008_tasmax"
+print(knmi_dataset)
+wrf311_dataset = local.load_file(FILE_2, "tasmax")
+wrf311_dataset.name = "AFRICA_UC-WRF311_CTL_ERAINT_MM_50km-rg_1989-2008_tasmax"
+print(wrf311_dataset)
+
+
+
+""" Step 2: Fetch an OCW Dataset Object from the data_source.rcmed module """
+print("Working with the rcmed interface to get CRU3.1 Daily-Max Temp")
+metadata = rcmed.get_parameters_metadata()
+
+cru_31 = [m for m in metadata if m['parameter_id'] == "39"][0]
+
+""" The RCMED API uses the following function to query, subset and return the 
+raw data from the database:
+
+rcmed.parameter_dataset(dataset_id, parameter_id, min_lat, max_lat, min_lon, 
+                        max_lon, start_time, end_time)
+
+The first two required params are in the cru_31 variable we defined earlier
+"""
+# Must cast to int since the rcmed api requires ints
+dataset_id = int(cru_31['dataset_id'])
+parameter_id = int(cru_31['parameter_id'])
+
+#  The spatial_boundaries() function returns the spatial extent of the dataset
+min_lat, max_lat, min_lon, max_lon = knmi_dataset.spatial_boundaries()
+
+print("Calculating the Maximum Overlap in Time for the datasets")
+
+cru_start = datetime.datetime.strptime(cru_31['start_date'], "%Y-%m-%d")
+cru_end = datetime.datetime.strptime(cru_31['end_date'], "%Y-%m-%d")
+knmi_start, knmi_end = knmi_dataset.time_range()
+# Grab the Max Start Time
+start_time = max([cru_start, knmi_start])
+# Grab the Min End Time
+end_time = min([cru_end, knmi_end])
+print("Overlap computed to be: %s to %s" % (start_time.strftime("%Y-%m-%d"),
+                                          end_time.strftime("%Y-%m-%d")))
+print("We are going to grab the first %s year(s) of data" % YEARS)
+end_time = datetime.datetime(start_time.year + YEARS, start_time.month, start_time.day)
+print("Final Overlap is: %s to %s" % (start_time.strftime("%Y-%m-%d"),
+                                          end_time.strftime("%Y-%m-%d")))
+
+
+print("Fetching data from RCMED...")
+cru31_dataset = rcmed.parameter_dataset(dataset_id,
+                                        parameter_id,
+                                        min_lat,
+                                        max_lat,
+                                        min_lon,
+                                        max_lon,
+                                        start_time,
+                                        end_time)
+
+import sys; sys.exit()
+""" Step 3: Resample Datasets so they are the same shape """
+print("CRU31_Dataset.values shape: (times, lats, lons) - %s" % (cru31_dataset.values.shape,))
+print("KNMI_Dataset.values shape: (times, lats, lons) - %s" % (knmi_dataset.values.shape,))
+print("Our two datasets have a mis-match in time. We will subset on time to %s years\n" %
YEARS)
+
+# Create a Bounds object to use for subsetting
+new_bounds = Bounds(min_lat, max_lat, min_lon, max_lon, start_time, end_time)
+knmi_dataset = dsp.subset(new_bounds, knmi_dataset)
+
+print("CRU31_Dataset.values shape: (times, lats, lons) - %s" % (cru31_dataset.values.shape,))
+print("KNMI_Dataset.values shape: (times, lats, lons) - %s \n" % (knmi_dataset.values.shape,))
+
+print("Temporally Rebinning the Datasets to a Single Timestep")
+# To run FULL temporal Rebinning use a timedelta > 366 days.  I used 999 in this example
+knmi_dataset = dsp.temporal_rebin(knmi_dataset, datetime.timedelta(days=999))
+cru31_dataset = dsp.temporal_rebin(cru31_dataset, datetime.timedelta(days=999))
+
+print("KNMI_Dataset.values shape: %s" % (knmi_dataset.values.shape,))
+print("CRU31_Dataset.values shape: %s \n\n" % (cru31_dataset.values.shape,))
+ 
+""" Spatially Regrid the Dataset Objects to a 1/2 degree grid """
+# Using the bounds we will create a new set of lats and lons on 1 degree step
+new_lons = np.arange(min_lon, max_lon, 0.5)
+new_lats = np.arange(min_lat, max_lat, 0.5)
+ 
+# Spatially regrid datasets using the new_lats, new_lons numpy arrays
+print("Spatially Regridding the KNMI_Dataset...")
+knmi_dataset = dsp.spatial_regrid(knmi_dataset, new_lats, new_lons)
+print("Spatially Regridding the CRU31_Dataset...")
+cru31_dataset = dsp.spatial_regrid(cru31_dataset, new_lats, new_lons)
+print("Final shape of the KNMI_Dataset:%s" % (knmi_dataset.values.shape, ))
+print("Final shape of the CRU31_Dataset:%s" % (cru31_dataset.values.shape, ))
+ 
+""" Step 4:  Build a Metric to use for Evaluation - Bias for this example """
+# You can build your own metrics, but OCW also ships with some common metrics
+print("Setting up a Bias metric to use for evaluation")
+bias = metrics.Bias()
+
+""" Step 5: Create an Evaluation Object using Datasets and our Metric """
+# The Evaluation Class Signature is:
+# Evaluation(reference, targets, metrics, subregions=None)
+# Evaluation can take in multiple targets and metrics, so we need to convert
+# our examples into Python lists.  Evaluation will iterate over the lists
+print("Making the Evaluation definition")
+bias_evaluation = evaluation.Evaluation(knmi_dataset, [cru31_dataset], [bias])
+print("Executing the Evaluation using the object's run() method")
+bias_evaluation.run()
+ 
+""" Step 6: Make a Plot from the Evaluation.results """
+# The Evaluation.results are a set of nested lists to support many different
+# possible Evaluation scenarios.
+#
+# The Evaluation results docs say:
+# The shape of results is (num_metrics, num_target_datasets) if no subregion
+# Accessing the actual results when we have used 1 metric and 1 dataset is
+# done this way:
+print("Accessing the Results of the Evaluation run")
+results = bias_evaluation.results[0][0]
+ 
+# From the bias output I want to make a Contour Map of the region
+print("Generating a contour map using ocw.plotter.draw_contour_map()")
+ 
+lats = new_lats
+lons = new_lons
+fname = OUTPUT_PLOT
+gridshape = (1, 1)  # Using a 1 x 1 since we have a single Bias for the full time range
+plot_title = "TASMAX Bias of KNMI Compared to CRU 3.1 (%s - %s)" % (start_time.strftime("%Y/%d/%m"),
end_time.strftime("%Y/%d/%m"))
+sub_titles = ["Full Temporal Range"]
+ 
+plotter.draw_contour_map(results, lats, lons, fname,
+                         gridshape=gridshape, ptitle=plot_title, 
+                         subtitles=sub_titles)


Mime
View raw message