USDA-ARS-NWRC / smrf

SMRF was designed to increase the flexibility of taking measured weather data, or atmospheric models, and distributing the data across a watershed.
Other
12 stars 4 forks source link

Local grid interpolation is slower than desired #91

Closed scotthavens closed 5 years ago

scotthavens commented 5 years ago

In the San Joaquin model run for a single day, we are seeing about 25% slow down due to the calculations in the local gradient method. Profiling the grid.py for the unit test test_load_data.test_grid_hrrr_local shows that about 90% of the time spent is in the calculation of the slopes and intercepts. The new interp_speed changes helps to speed up the interpolation but the test case is too small to truly see the speed increases.

micah-prime commented 5 years ago

This is the result of the line profiler over one timestep in the Tuolumne, which should be a good representation of where the code is spending it's time. This was using the cubic method for interpolation, which, unlike the linear method, does not show a speedup when changing from griddata to grid_interpolate_deconstructed.

It's evident here that the biggest offenders are:

122       852    1121563.0   1316.4     13.4              elev = pv.loc[drow].elevation
123       852    1395132.0   1637.5     16.7              pp = np.polyfit(elev, pv.loc[drow].data, 1)
130       852    1807012.0   2120.9     21.6              pv.loc[grid_name, ['slope', 'intercept']] = pp

138         1     635758.0 635758.0      7.6          grid_slope_og = griddata((pv.utm_x, pv.utm_y), pv.slope, (self.GridX, self.GridY), method=grid_method)
139         1     615866.0 615866.0      7.4          grid_slope = grid_interpolate_deconstructed(self.tri, pv.slope.values[:], (self.GridX, self.GridY), method=grid_method)
142         1     654457.0 654457.0      7.8          grid_intercept_og = griddata((pv.utm_x, pv.utm_y), pv.intercept, (self.GridX, self.GridY), method=grid_method)
143         1     626213.0 626213.0      7.5          grid_intercept = grid_interpolate_deconstructed(self.tri, pv.intercept.values[:], (self.GridX, self.GridY), method=grid_method)
151         1     655289.0 655289.0      7.8          idtrend_og = griddata((pv.utm_x, pv.utm_y), dtrend, (self.GridX, self.GridY), method=grid_method)
152         1     620712.0 620712.0      7.4          idtrend = grid_interpolate_deconstructed(self.tri, dtrend, (self.GridX, self.GridY), method=grid_method)

The griddata functions account for ~22% of the time. The pandas data grabs and the polyfit are quite slow.

Here is the full print out:

Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
Timer unit: 1e-06 s

Total time: 8.34768 s
File: /home/micahsandusky/Documents/Code/virtualenvirons/env3/lib/python3.5/site-packages/smrf-0.8.8-py3.5-linux-x86_64.egg/smrf/spatial/grid.py
Function: detrendedInterpolationLocal at line 105

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   105                                               def detrendedInterpolationLocal(self, data, flag=0, grid_method='linear'):
   106                                                   """
   107                                                   Interpolate using a detrended approach
   108                                           
   109                                                   Args:
   110                                                       data: data to interpolate
   111                                                       grid_method: scipy.interpolate.griddata interpolation method
   112                                                   """
   113                                           
   114                                                   # copy the metadata to append the slope/intercept and data to
   115         1        254.0    254.0      0.0          pv = self.metadata.copy()
   116         1        972.0    972.0      0.0          pv['slope'] = np.nan
   117         1        903.0    903.0      0.0          pv['intercept'] = np.nan
   118         1        919.0    919.0      0.0          pv['data'] = data
   119                                           
   120       853     191236.0    224.2      2.3          for grid_name,drow in self.dist_df.iterrows():
   121                                           
   122       852    1121563.0   1316.4     13.4              elev = pv.loc[drow].elevation
   123       852    1395132.0   1637.5     16.7              pp = np.polyfit(elev, pv.loc[drow].data, 1)
   124                                           
   125                                                       # apply trend constraints
   126       852       2690.0      3.2      0.0              if flag == 1 and pp[0] < 0:
   127       186        968.0      5.2      0.0                  pp = np.array([0, 0])
   128       666        602.0      0.9      0.0              elif (flag == -1 and pp[0] > 0):
   129                                                           pp = np.array([0, 0])
   130       852    1807012.0   2120.9     21.6              pv.loc[grid_name, ['slope', 'intercept']] = pp
   131                                           
   132                                                   # get triangulation
   133         1          1.0      1.0      0.0          if self.tri is None:
   134                                                       xy = _ndim_coords_from_arrays((pv.utm_x, pv.utm_y))
   135                                                       self.tri = qhull.Delaunay(xy)
   136                                           
   137                                                   # interpolate the slope/intercept
   138         1     635758.0 635758.0      7.6          grid_slope_og = griddata((pv.utm_x, pv.utm_y), pv.slope, (self.GridX, self.GridY), method=grid_method)
   139         1     615866.0 615866.0      7.4          grid_slope = grid_interpolate_deconstructed(self.tri, pv.slope.values[:], (self.GridX, self.GridY), method=grid_method)
   140         1       2981.0   2981.0      0.0          print('Slope difference: {}'.format(np.sum(grid_slope_og-grid_slope)))
   141                                           
   142         1     654457.0 654457.0      7.8          grid_intercept_og = griddata((pv.utm_x, pv.utm_y), pv.intercept, (self.GridX, self.GridY), method=grid_method)
   143         1     626213.0 626213.0      7.5          grid_intercept = grid_interpolate_deconstructed(self.tri, pv.intercept.values[:], (self.GridX, self.GridY), method=grid_method)
   144         1       3084.0   3084.0      0.0          print('Intercept difference: {}'.format(np.sum(grid_intercept_og-grid_intercept)))
   145                                           
   146                                                   # remove the elevation trend from the HRRR precip
   147         1        788.0    788.0      0.0          el_trend = pv.elevation * pv.slope + pv.intercept
   148         1        361.0    361.0      0.0          dtrend = pv.data - el_trend
   149                                           
   150                                                   # interpolate the residuals over the DEM
   151         1     655289.0 655289.0      7.8          idtrend_og = griddata((pv.utm_x, pv.utm_y), dtrend, (self.GridX, self.GridY), method=grid_method)
   152         1     620712.0 620712.0      7.4          idtrend = grid_interpolate_deconstructed(self.tri, dtrend, (self.GridX, self.GridY), method=grid_method)
   153         1       3563.0   3563.0      0.0          print('idtrend difference: {}'.format(np.sum(idtrend_og-idtrend)))
   154                                           
   155                                                   # reinterpolate
   156         1       6350.0   6350.0      0.1          rtrend = idtrend + grid_slope * self.GridZ + grid_intercept
   157                                           
   158         1          3.0      3.0      0.0          return rtrend
micah-prime commented 5 years ago

Here is the timing with the comparison of the Pandas changes

Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
Timer unit: 1e-06 s

Total time: 8.26371 s
File: /home/micahsandusky/Documents/Code/virtualenvirons/env3/lib/python3.5/site-packages/smrf-0.8.8-py3.5-linux-x86_64.egg/smrf/spatial/grid.py
Function: detrendedInterpolationLocal at line 118

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   118                                               def detrendedInterpolationLocal(self, data, flag=0, grid_method='linear'):
   119                                                   """
   120                                                   Interpolate using a detrended approach
   121                                           
   122                                                   Args:
   123                                                       data: data to interpolate
   124                                                       grid_method: scipy.interpolate.griddata interpolation method
   125                                                   """
   126                                           
   127                                                   # copy the metadata to append the slope/intercept and data to
   128         1        390.0    390.0      0.0          pv = self.metadata.copy()
   129         1        815.0    815.0      0.0          pv['slope'] = np.nan
   130         1        660.0    660.0      0.0          pv['intercept'] = np.nan
   131         1        630.0    630.0      0.0          pv['data'] = data
   132                                           
   133                                                   # take the new full_df and fill a data column
   134                                                   # Adapted from:
   135                                                   # https://stackoverflow.com/questions/51140302/using-polyfit-on-pandas-dataframe-and-then-adding-the-results-to-new-columns
   136         1        552.0    552.0      0.0          df = self.full_df.copy()
   137         1       5694.0   5694.0      0.1          df['data'] = data[df['cell_local']].values
   138         1       2199.0   2199.0      0.0          df = df.set_index('cell_id')
   139         1     318227.0 318227.0      3.9          df['fit'] = df.groupby('cell_id').apply(lambda x: np.polyfit(x.elevation, x.data, 1))
   140                                           
   141                                                   # drop all the duplicates
   142         1        768.0    768.0      0.0          df.reset_index(inplace=True)
   143         1       2861.0   2861.0      0.0          df.drop_duplicates(subset=['cell_id'], keep='first', inplace=True)
   144         1        699.0    699.0      0.0          df.set_index('cell_id', inplace=True)
   145         1     157529.0 157529.0      1.9          df[['slope', 'intercept']] = df.fit.apply(pd.Series)
   146                                                   # df = df.drop(columns='fit').reset_index()
   147                                           
   148       853     124013.0    145.4      1.5          for grid_name,drow in self.dist_df.iterrows():
   149                                           
   150       852    1056591.0   1240.1     12.8              elev = pv.loc[drow].elevation
   151       852    1242339.0   1458.1     15.0              pp = np.polyfit(elev, pv.loc[drow].data, 1)
   152                                           
   153                                                       # apply trend constraints
   154       852       1518.0      1.8      0.0              if flag == 1 and pp[0] < 0:
   155                                                           pp = np.array([0, 0])
   156       852       2367.0      2.8      0.0              elif (flag == -1 and pp[0] > 0):
   157                                                           pp = np.array([0, 0])
   158       852    1149725.0   1349.4     13.9              pv.loc[grid_name, ['slope', 'intercept']] = pp
   159                                           
   160                                                   # compare the pv and df
   161         1        667.0    667.0      0.0          print('Cell slope difference: {}'.format(np.mean(df['slope']-pv['slope'])))
   162         1        625.0    625.0      0.0          print('Cell intercept difference: {}'.format(np.mean(df['intercept']-pv['intercept'])))
   163                                           
   164                                                   # get triangulation
   165         1          6.0      6.0      0.0          if self.tri is None:
   166                                                       xy = _ndim_coords_from_arrays((pv.utm_x, pv.utm_y))
   167                                                       self.tri = qhull.Delaunay(xy)
   168                                           
   169                                                   # interpolate the slope/intercept
   170         1     716195.0 716195.0      8.7          grid_slope_og = griddata((pv.utm_x, pv.utm_y), pv.slope, (self.GridX, self.GridY), method=grid_method)
   171         1     615731.0 615731.0      7.5          grid_slope = grid_interpolate_deconstructed(self.tri, pv.slope.values[:], (self.GridX, self.GridY), method=grid_method)
   172         1       5433.0   5433.0      0.1          print('Slope difference: {}'.format(np.sum(grid_slope_og-grid_slope)))
   173                                           
   174         1     638040.0 638040.0      7.7          grid_intercept_og = griddata((pv.utm_x, pv.utm_y), pv.intercept, (self.GridX, self.GridY), method=grid_method)
   175         1     761714.0 761714.0      9.2          grid_intercept = grid_interpolate_deconstructed(self.tri, pv.intercept.values[:], (self.GridX, self.GridY), method=grid_method)
   176         1       6076.0   6076.0      0.1          print('Intercept difference: {}'.format(np.sum(grid_intercept_og-grid_intercept)))
   177                                           
   178                                                   # remove the elevation trend from the HRRR precip
   179         1        669.0    669.0      0.0          el_trend = pv.elevation * pv.slope + pv.intercept
   180         1        287.0    287.0      0.0          dtrend = pv.data - el_trend
   181                                           
   182                                                   # interpolate the residuals over the DEM
   183         1     717946.0 717946.0      8.7          idtrend_og = griddata((pv.utm_x, pv.utm_y), dtrend, (self.GridX, self.GridY), method=grid_method)
   184         1     724089.0 724089.0      8.8          idtrend = grid_interpolate_deconstructed(self.tri, dtrend, (self.GridX, self.GridY), method=grid_method)
   185         1       3065.0   3065.0      0.0          print('idtrend difference: {}'.format(np.sum(idtrend_og-idtrend)))
   186                                           
   187                                                   # reinterpolate
   188         1       5586.0   5586.0      0.1          rtrend = idtrend + grid_slope * self.GridZ + grid_intercept
   189                                           
   190         1          4.0      4.0      0.0          return rtrend
scotthavens commented 5 years ago

So with the new vectorized pandas instead of the for loop, the speed up appears to be from about 40% of the time to now 6% of the time when all are done in the same loop.

Through the Tuolumne test, have we confirmed that there are no differences between the two methods?

micah-prime commented 5 years ago

I think the print out at the beginning of the last timing log shows that, I am correct there?

micah-prime commented 5 years ago

On second thought, it does appear to be a little different on some timesteps:

DEBUG:smrf.distribute.air_temp:2019-04-02 05:00:00+00:00 -- Distributing air_temp
Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
DEBUG:smrf.distribute.vapor_pressure:2019-04-02 05:00:00+00:00 -- Distributing vapor_pressure
Cell slope difference: 5.123867052388005e-07
Cell intercept difference: 0.563147768546374
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
DEBUG:smrf.distribute.vapor_pressure:2019-04-02 05:00:00+00:00 -- Calculating dew point
DEBUG:smrf.distribute.wind:2019-04-02 05:00:00+00:00 Distributing wind_direction and wind_speed
DEBUG:smrf.distribute.precipitation:2019-04-02 05:00:00+00:00 Distributing all precip
DEBUG:smrf.distribute.precipitation:Storming? True
DEBUG:smrf.distribute.precipitation:Current Storm ID = 0
Cell slope difference: -1.5027324964643772e-05
Cell intercept difference: 0.08608602142206176
Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
micah-prime commented 5 years ago

Looking good with the latest fix! (constraining the fit)

DEBUG:smrf.distribute.air_temp:2019-04-02 05:00:00+00:00 -- Distributing air_temp
Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
DEBUG:smrf.distribute.vapor_pressure:2019-04-02 05:00:00+00:00 -- Distributing vapor_pressure
Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
DEBUG:smrf.distribute.vapor_pressure:2019-04-02 05:00:00+00:00 -- Calculating dew point
DEBUG:smrf.distribute.wind:2019-04-02 05:00:00+00:00 Distributing wind_direction and wind_speed
DEBUG:smrf.distribute.precipitation:2019-04-02 05:00:00+00:00 Distributing all precip
DEBUG:smrf.distribute.precipitation:Storming? True
DEBUG:smrf.distribute.precipitation:Current Storm ID = 0
Cell slope difference: 0.0
Cell intercept difference: 0.0
Slope difference: 0.0
Intercept difference: 0.0
idtrend difference: 0.0
micah-prime commented 5 years ago

It looks like the time per timestep is way down and the gridded interpolation is by far the slowest part, which is how it should be. This issue is just about remedied.

Timer unit: 1e-06 s

Total time: 2.62898 s
File: /home/micahsandusky/Documents/Code/virtualenvirons/env3/lib/python3.5/site-packages/smrf-0.8.8-py3.5-linux-x86_64.egg/smrf/spatial/grid.py
Function: detrendedInterpolationLocal at line 118

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   118                                               def detrendedInterpolationLocal(self, data, flag=0, grid_method='linear'):
   119                                                   """
   120                                                   Interpolate using a detrended approach
   121                                           
   122                                                   Args:
   123                                                       data: data to interpolate
   124                                                       grid_method: scipy.interpolate.griddata interpolation method
   125                                                   """
   126                                           
   127                                                   # take the new full_df and fill a data column
   128                                                   # Adapted from:
   129                                                   # https://stackoverflow.com/questions/51140302/using-polyfit-on-pandas-dataframe-and-then-adding-the-results-to-new-columns
   130         1        552.0    552.0      0.0          df = self.full_df.copy()
   131         1       4376.0   4376.0      0.2          df['data'] = data[df['cell_local']].values
   132         1       1951.0   1951.0      0.1          df = df.set_index('cell_id')
   133         1     385707.0 385707.0     14.7          df['fit'] = df.groupby('cell_id').apply(lambda x: np.polyfit(x.elevation, x.data, 1))
   134                                           
   135                                                   # drop all the duplicates
   136         1       1195.0   1195.0      0.0          df.reset_index(inplace=True)
   137         1       3714.0   3714.0      0.1          df.drop_duplicates(subset=['cell_id'], keep='first', inplace=True)
   138         1        814.0    814.0      0.0          df.set_index('cell_id', inplace=True)
   139         1     205773.0 205773.0      7.8          df[['slope', 'intercept']] = df.fit.apply(pd.Series)
   140                                                   # df = df.drop(columns='fit').reset_index()
   141                                           
   142                                                   # apply trend constraints
   143         1          2.0      2.0      0.0          if flag == 1:
   144                                                       df.loc[df['slope'] < 0, ['slope', 'intercept']] = 0
   145         1          1.0      1.0      0.0          elif flag == -1:
   146         1       7133.0   7133.0      0.3              df.loc[df['slope'] > 0, ['slope', 'intercept']] = 0
   147                                           
   148                                                   # get triangulation
   149         1          6.0      6.0      0.0          if self.tri is None:
   150                                                       xy = _ndim_coords_from_arrays((self.metadata.utm_x, self.metadata.utm_y))
   151                                                       self.tri = qhull.Delaunay(xy)
   152                                           
   153                                                   # interpolate the slope/intercept
   154         1     690995.0 690995.0     26.3          grid_slope = grid_interpolate_deconstructed(self.tri, df.slope.values[:], (self.GridX, self.GridY), method=grid_method)
   155                                           
   156         1     642615.0 642615.0     24.4          grid_intercept = grid_interpolate_deconstructed(self.tri, df.intercept.values[:], (self.GridX, self.GridY), method=grid_method)
   157                                           
   158                                                   # remove the elevation trend from the HRRR precip
   159         1        684.0    684.0      0.0          el_trend = df.elevation * df.slope + df.intercept
   160         1        272.0    272.0      0.0          dtrend = df.data - el_trend
   161                                           
   162                                                   # interpolate the residuals over the DEM
   163         1     676623.0 676623.0     25.7          idtrend = grid_interpolate_deconstructed(self.tri, dtrend, (self.GridX, self.GridY), method=grid_method)
   164                                           
   165                                                   # reinterpolate
   166         1       6559.0   6559.0      0.2          rtrend = idtrend + grid_slope * self.GridZ + grid_intercept
   167                                           
   168         1          3.0      3.0      0.0          return rtrend