Skip to content
This repository has been archived by the owner. It is now read-only.
Permalink
Browse files
Merge branch 'CLIMATE-841' of https://github.com/huikyole/climate
  • Loading branch information
agoodm committed Aug 10, 2016
2 parents 625eec2 + 243a895 commit 5d1bc5c03a70ee196ef0011ee8e22fc6c9afe48d
Show file tree
Hide file tree
Showing 5 changed files with 182 additions and 82 deletions.
File renamed without changes.
@@ -0,0 +1,56 @@
# 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.

#Apache OCW lib immports
import ocw.dataset_processor as dsp
import ocw.utils as utils
from ocw.dataset import Bounds
import ocw.data_source.rcmed as rcmed
import ocw.plotter as plotter

from datetime import datetime
import numpy.ma as ma

import ssl

if hasattr(ssl, '_create_unverified_context'):
ssl._create_default_https_context = ssl._create_unverified_context

# rectangular boundary
min_lat = 15.75
max_lat = 55.75
min_lon = -125.75
max_lon = -66.75

start_time = datetime(1998,1,1)
end_time = datetime(1998,12,31)

TRMM_dataset = rcmed.parameter_dataset(3, 36, min_lat, max_lat, min_lon, max_lon,
start_time, end_time)

Cuba_and_Bahamas_bounds = Bounds(boundary_type='countries', countries=['Cuba','Bahamas'])
TRMM_dataset2 = dsp.subset(TRMM_dataset, Cuba_and_Bahamas_bounds, extract=False) # to mask out the data over Mexico and Canada

plotter.draw_contour_map(ma.mean(TRMM_dataset2.values, axis=0), TRMM_dataset2.lats, TRMM_dataset2.lons, fname='TRMM_without_Cuba_and_Bahamas')

NCA_SW_bounds = Bounds(boundary_type='us_states', us_states=['CA','NV','UT','AZ','NM','CO'])
TRMM_dataset3 = dsp.subset(TRMM_dataset2, NCA_SW_bounds, extract=True) # to mask out the data over Mexico and Canada

plotter.draw_contour_map(ma.mean(TRMM_dataset3.values, axis=0), TRMM_dataset3.lats, TRMM_dataset3.lons, fname='TRMM_NCA_SW')



@@ -298,9 +298,9 @@ def __init__(self, boundary_type='rectangular',
self._end = None

if boundary_type == 'us_states':
self.masked_regions = shapefile_boundary(boundary_type, us_states)
self.masked_regions = utils.shapefile_boundary(boundary_type, us_states)
if boundary_type == 'countries':
self.masked_regions = shapefile_boundary(boundary_type, countries)
self.masked_regions = utils.shapefile_boundary(boundary_type, countries)
if boundary_type == 'user':
file_object = netCDF4.Dataset(user_mask_file)
self.mask_variable = file_object.variables[mask_variable_name][:]
@@ -334,7 +334,7 @@ def __init__(self, boundary_type='rectangular',
self.lon_min = float(lon_min)
self.lon_max = float(lon_max)
if boundary_type[:6].upper() == 'CORDEX':
self.lat_min, self.lat_max, self.lon_min, self.lon_max = CORDEX_boundary(boundary_type[6:].replace(" ","").lower())
self.lat_min, self.lat_max, self.lon_min, self.lon_max = utils.CORDEX_boundary(boundary_type[6:].replace(" ","").lower())

@property
def start(self):
@@ -364,62 +364,3 @@ def end(self, value):

self._end = value

def shapefile_boundary(boundary_type, region_names):
'''
:param boundary_type: The type of spatial subset boundary
:type boundary_type: :mod:'string'
:param region_names: An array of regions for spatial subset
:type region_names: :mod:'list'
'''
# Read the shapefile
map_read = Basemap()
regions = []
shapefile_dir = os.sep.join([os.path.dirname(__file__), 'shape'])
map_read.readshapefile(os.path.join(shapefile_dir, boundary_type),
boundary_type)
if boundary_type == 'us_states':
for region_name in region_names:
for iregion, region_info in enumerate(map_read.us_states_info):
if region_info['st'] == region_name:
regions.append(numpy.array(map_read.us_states[iregion]))
elif boundary_type == 'countries':
for region_name in region_names:
for iregion, region_info in enumerate(map_read.countries_info):
if region_info['COUNTRY'] == region_name:
regions.append(numpy.array(map_read.countries[iregion]))
return regions

def CORDEX_boundary(domain_name):
'''
:param domain_name: CORDEX domain name (http://www.cordex.org/)
:type domain_name: :mod:'string'
'''
if domain_name =='southamerica':
return -57.61, 18.50, 254.28-360., 343.02-360.
if domain_name =='centralamerica':
return -19.46, 34.83, 235.74-360., 337.78-360.
if domain_name =='northamerica':
return 12.55, 75.88, 189.26-360., 336.74-360.
if domain_name =='europe':
return 22.20, 71.84, 338.23-360., 64.4
if domain_name =='africa':
return -45.76, 42.24, 335.36-360., 60.28
if domain_name =='southasia':
return -15.23, 45.07, 19.88, 115.55
if domain_name =='eastasia':
return -0.10, 61.90, 51.59, 179.99
if domain_name =='centralasia':
return 18.34, 69.37, 11.05, 139.13
if domain_name =='australasia':
return -52.36, 12.21, 89.25, 179.99
if domain_name =='antartica':
return -89.48,-56.00, -179.00, 179.00
if domain_name =='artic':
return 46.06, 89.50, -179.00, 179.00
if domain_name =='mediterranean':
return 25.63, 56.66, 339.79-360.00, 50.85
if domain_name =='middleeastnorthafrica':
return -7.00, 45.00, 333.00-360.00, 76.00
if domain_name =='southeastasia':
return -15.14, 27.26, 89.26, 146.96
@@ -16,6 +16,7 @@
#

from ocw import dataset as ds
import ocw.utils as utils

import datetime
import numpy as np
@@ -362,7 +363,7 @@ def ensemble(datasets):
return ensemble_dataset


def subset(target_dataset, subregion, subregion_name=None):
def subset(target_dataset, subregion, subregion_name=None, extract=True):
'''Subset given dataset(s) with subregion information
:param subregion: The Bounds with which to subset the target Dataset.
@@ -374,6 +375,9 @@ def subset(target_dataset, subregion, subregion_name=None):
:param subregion_name: The subset-ed Dataset name
:type subregion_name: :mod:`string`
:param extract: If False, the dataset inside regions will be masked.
:type extract: :mod:`boolean`
:returns: The subset-ed Dataset object
:rtype: :class:`dataset.Dataset`
@@ -387,7 +391,7 @@ def subset(target_dataset, subregion, subregion_name=None):
if not subregion_name:
subregion_name = target_dataset.name

if subregion.lat_min:
if hasattr(subregion, 'lat_min'):
# If boundary_type is 'rectangular' or 'CORDEX ***', ensure that the subregion information is well formed
_are_bounds_contained_by_dataset(target_dataset, subregion)

@@ -435,27 +439,36 @@ def subset(target_dataset, subregion, subregion_name=None):
dataset_slices["time_start"]:dataset_slices["time_end"] + 1,
dataset_slices["lat_start"]:dataset_slices["lat_end"] + 1,
dataset_slices["lon_start"]:dataset_slices["lon_end"] + 1]

# Build new dataset with subset information
return ds.Dataset(
# Build new dataset with subset information
return ds.Dataset(
# Slice the lats array with our calculated slice indices
target_dataset.lats[dataset_slices["lat_start"]:
dataset_slices["lat_end"] + 1],
target_dataset.lats[dataset_slices["lat_start"]:
dataset_slices["lat_end"] + 1],
# Slice the lons array with our calculated slice indices
target_dataset.lons[dataset_slices["lon_start"]:
dataset_slices["lon_end"] + 1],
target_dataset.lons[dataset_slices["lon_start"]:
dataset_slices["lon_end"] + 1],
# Slice the times array with our calculated slice indices
target_dataset.times[dataset_slices["time_start"]:
dataset_slices["time_end"] + 1],
target_dataset.times[dataset_slices["time_start"]:
dataset_slices["time_end"] + 1],
# Slice the values array with our calculated slice indices
subset_values,
variable=target_dataset.variable,
units=target_dataset.units,
name=subregion_name,
origin=target_dataset.origin
)


subset_values,
variable=target_dataset.variable,
units=target_dataset.units,
name=subregion_name,
origin=target_dataset.origin
)

if subregion.boundary_type == 'us_states' or subregion.boundary_type == 'countries':
start_time_index = np.where(
target_dataset.times == subregion.start)[0][0]
end_time_index = np.where(target_dataset.times == subregion.end)[0][0]
target_dataset = temporal_slice(
target_dataset, start_time_index, end_time_index)
nt, ny, nx = target_dataset.values.shape
spatial_mask = utils.mask_using_shapefile_info(target_dataset.lons, target_dataset.lats, subregion.masked_regions, extract = extract)
target_dataset.values = utils.propagate_spatial_mask_over_time(target_dataset.values, mask=spatial_mask)
return target_dataset

def temporal_slice(target_dataset, start_time_index, end_time_index):
'''Temporally slice given dataset(s) with subregion information. This does not
spatially subset the target_Dataset
@@ -18,11 +18,13 @@
''''''

import sys
import os
import datetime as dt
import numpy as np
import numpy.ma as ma

from mpl_toolkits.basemap import shiftgrid
from mpl_toolkits.basemap import shiftgrid, Basemap, maskoceans
from matplotlib.path import Path
from dateutil.relativedelta import relativedelta
from netCDF4 import num2date

@@ -445,3 +447,91 @@ def calc_area_weighted_spatial_average(dataset, area_weight=False):
spatial_average[it] = ma.average(dataset.values[it, :])

return spatial_average

def shapefile_boundary(boundary_type, region_names):
'''
:param boundary_type: The type of spatial subset boundary
:type boundary_type: :mod:'string'
:param region_names: An array of regions for spatial subset
:type region_names: :mod:'list'
'''
# Read the shapefile
map_read = Basemap()
regions = []
shapefile_dir = os.sep.join([os.path.dirname(__file__), 'shape'])
map_read.readshapefile(os.path.join(shapefile_dir, boundary_type),
boundary_type)
if boundary_type == 'us_states':
for region_name in region_names:
for iregion, region_info in enumerate(map_read.us_states_info):
if region_info['st'] == region_name:
regions.append(np.array(map_read.us_states[iregion]))
elif boundary_type == 'countries':
for region_name in region_names:
for iregion, region_info in enumerate(map_read.countries_info):
if region_info['COUNTRY'].replace(" ","").lower() == region_name.replace(" ","").lower():
regions.append(np.array(map_read.countries[iregion]))
return regions

def CORDEX_boundary(domain_name):
'''
:param domain_name: CORDEX domain name (http://www.cordex.org/)
:type domain_name: :mod:'string'
'''
if domain_name =='southamerica':
return -57.61, 18.50, 254.28-360., 343.02-360.
if domain_name =='centralamerica':
return -19.46, 34.83, 235.74-360., 337.78-360.
if domain_name =='northamerica':
return 12.55, 75.88, 189.26-360., 336.74-360.
if domain_name =='europe':
return 22.20, 71.84, 338.23-360., 64.4
if domain_name =='africa':
return -45.76, 42.24, 335.36-360., 60.28
if domain_name =='southasia':
return -15.23, 45.07, 19.88, 115.55
if domain_name =='eastasia':
return -0.10, 61.90, 51.59, 179.99
if domain_name =='centralasia':
return 18.34, 69.37, 11.05, 139.13
if domain_name =='australasia':
return -52.36, 12.21, 89.25, 179.99
if domain_name =='antartica':
return -89.48,-56.00, -179.00, 179.00
if domain_name =='artic':
return 46.06, 89.50, -179.00, 179.00
if domain_name =='mediterranean':
return 25.63, 56.66, 339.79-360.00, 50.85
if domain_name =='middleeastnorthafrica':
return -7.00, 45.00, 333.00-360.00, 76.00
if domain_name =='southeastasia':
return -15.14, 27.26, 89.26, 146.96

def mask_using_shapefile_info(lons, lats, masked_regions, extract = True):
if lons.ndim == 2 and lats.ndim == 2:
lons_2d = lons
lats_2d = lats
else:
lons_2d, lats_2d = np.meshgrid(lons, lats)

points = np.array((lons_2d.flatten(), lats_2d.flatten())).T
for iregion, region in enumerate(masked_regions):
mpath = Path(region)
mask0 = mpath.contains_points(points).reshape(lons_2d.shape)
if iregion == 0:
mask = mask0
else:
mask = mask|mask0
if extract:
mask = np.invert(mask)
return mask

def propagate_spatial_mask_over_time(data_array, mask):
if data_array.ndim == 3 and mask.ndim == 2:
nt = data_array.shape[0]
for it in np.arange(nt):
mask_original = data_array[it,:].mask
data_array[it,:] = ma.array(data_array[it,:], mask=mask|mask_original)
return data_array

0 comments on commit 5d1bc5c

Please sign in to comment.