Index: trunk/SDToolBox/output_messages.py =================================================================== diff -u -r59 -r60 --- trunk/SDToolBox/output_messages.py (.../output_messages.py) (revision 59) +++ trunk/SDToolBox/output_messages.py (.../output_messages.py) (revision 60) @@ -21,3 +21,5 @@ error_all_arguments_required = 'All arguments are required for resampling.' error_max_lon_not_set = 'Max longitud needs to be set.' error_needs_to_be_in_subclass = 'Needs to be implemented in subclasses.' +error_no_gradient_was_calculated = 'No gradient was calculated, ' + \ + ' please calculate the gradient on the xarray before proceeding' Index: trunk/tests/test_data_processing.py =================================================================== diff -u -r58 -r60 --- trunk/tests/test_data_processing.py (.../test_data_processing.py) (revision 58) +++ trunk/tests/test_data_processing.py (.../test_data_processing.py) (revision 60) @@ -138,3 +138,50 @@ # 3. Then assert str(e_info.value) == expected_message + + +class TestSpatialGradientsCalculation: + + @pytest.mark.unittest + def test_given_an_xarray_compute_slp_gradients(self): + lat = [43.125, 50, 60.125, 13.5] + lon = [43.125, 50, 60.125, 13.5] + times = pd.date_range('2000-01-01', periods=4) + coordinates = ['lon', 'lat'] + + data = np.random.rand(len(times), len(lon), len(lat)) + + data_array = xr.DataArray( + data, + coords=[times, lon, lat], + dims=['times', 'longitude', 'latitude']) + assert data_array is not None + + data_processing = DataProcessing() + result = data_processing.compute_spatial_gradients(data_array) + + +class TestComputePCA: + + @pytest.mark.unittest + def test_given_an_xarray_compute_PCA(self): + lat = [43.125, 50, 60.125] + lon = [43.125, 50, 60.125] + times = pd.date_range('2000-01-01', periods=3) + coordinates = ['lon', 'lat'] + + data = np.random.rand(len(times), len(lon), len(lat)) + + data_array = xr.DataArray( + data, + coords=[times, lon, lat], + dims=['times', 'longitude', 'latitude']) + assert data_array is not None + + # data_processing = DataProcessing() + # result = data_processing.compute_spatial_gradients(data_array) + + # data = [] + # grads = xr.DataArray( + + # ) Index: trunk/SDToolBox/data_processing.py =================================================================== diff -u -r58 -r60 --- trunk/SDToolBox/data_processing.py (.../data_processing.py) (revision 58) +++ trunk/SDToolBox/data_processing.py (.../data_processing.py) (revision 60) @@ -11,6 +11,8 @@ from SDToolBox import output_messages as om from SDToolBox.output_data import OutputData +from sklearn.decomposition import PCA as skPCA +from sklearn.preprocessing import StandardScaler import numpy as np import xarray as xr @@ -87,6 +89,72 @@ frequency_string=frequency_string ) - @staticmethod - def spatial_resampling(): - pass + def compute_spatial_gradients(self, slp: xr.DataArray): + + meshed_latitudes, meshed_longitudes = \ + np.meshgrid(slp.coords['latitude'], slp.coords['longitude']) + + times = slp.coords['times'] + latitude = slp.coords['latitude'] + longitude = slp.coords['longitude'] + + latitude_size = len(latitude) + longitude_size = len(longitude) + times_size = len(times) + + derivative_points = \ + np.zeros((times_size, latitude_size, longitude_size), dtype=float) + + gradient_slp = [] + + for time in range(1, len(times)): + for longitude_idx in range(2, len(longitude)-1): # x + for latitude_idx in range(2, len(latitude)-1): # y + + phi = np.pi * \ + np.abs(meshed_latitudes[latitude_idx, longitude_idx])/180 + derivative_longitude_first = \ + (slp[:, latitude_idx, longitude_idx] - slp[:, latitude_idx, longitude_idx - 1]) / np.cos(phi) + derivative_longitude_second = \ + (slp[:, latitude_idx, longitude_idx + 1] - slp[:, latitude_idx, longitude_idx]) / np.cos(phi) + derivative_latitude_first = \ + (slp[:, latitude_idx, longitude_idx] - slp[:, latitude_idx - 1, longitude_idx]) / np.cos(phi) + derivative_latitude_second = \ + (slp[:, latitude_idx + 1, longitude_idx] - slp[:, latitude_idx, longitude_idx]) / np.cos(phi) + + derivative_points[:, latitude_idx, longitude_idx] = np.square(derivative_longitude_first) + np.square(derivative_longitude_second)/2 + np.square(derivative_latitude_first) + np.square(derivative_latitude_second)/2 + gradient_slp.append(derivative_points) + # probably here we need to use .assign_coords function to add grad2slp to the existing array + + data_array = xr.DataArray( + gradient_slp, + coords=[time, latitude, longitude], + dims=['times', 'longitude', 'latitude']) + + return data_array + + def compute_PCA(self, slp: xr.DataArray): + + grad2slp = slp.coords['grad2slp'] + + if grad2slp is None: + raise Exception(om.error_no_gradient_was_calculated) + + scaler = StandardScaler(with_mean=True, with_std=True) + scaler.fit(grad2slp) + transformed_grad2slp = scaler.transform(grad2slp) + exvar = 99.88/100 + pca = skPCA(n_components=exvar) + pca.fit(transformed_grad2slp) + principal_components = pca.components_ + temporal_indices = \ + np.dot(transformed_grad2slp, principal_components.transpose()) + standard_temporal_indices = temporal_indices.std(axis=0) + pca.mean = scaler.mean_() + pca.std = np.std(scaler.var_()) + pca.eof = principal_components + pca.pc = temporal_indices + pca.variance = \ + standard_temporal_indices**2/np.sum(standard_temporal_indices**2)*100 + pca.variance = \ + pca.explained_variance_ratio_() Index: trunk/SDToolBox/extract_data_EARTH.py =================================================================== diff -u -r58 -r60 --- trunk/SDToolBox/extract_data_EARTH.py (.../extract_data_EARTH.py) (revision 58) +++ trunk/SDToolBox/extract_data_EARTH.py (.../extract_data_EARTH.py) (revision 60) @@ -47,9 +47,13 @@ for n_variable, variable_name in enumerate(filtered_dict): for year in self._input_years: if month < 10: - base_file_name = self._input_EARTH_scenario + '_\\' + 'EC-Earth_RCP4.5_MSLP_' + str(year) + '0' + str(month) + '.nc' + base_file_name = self._input_EARTH_scenario + \ + '_\\' + 'EC-Earth_RCP4.5_MSLP_' + \ + str(year) + '0' + str(month) + '.nc' else: - base_file_name = self._input_EARTH_scenario + '_\\' + 'EC-Earth_RCP4.5_MSLP_' + str(year) + str(month) + '.nc' + base_file_name = self._input_EARTH_scenario + \ + '_\\' + 'EC-Earth_RCP4.5_MSLP_' + \ + str(year) + str(month) + '.nc' # If file does not exist simply go to the next one # case_dir = os.path.join(directory_path, case_name_value) case_file_path = os.path.join(directory_path, base_file_name)