HARPgroup / model_meteorology

0 stars 0 forks source link

Reversibility of Spatial Resample and Temporal Disaggregate Steps #79

Open rburghol opened 1 month ago

rburghol commented 1 month ago

It appears that an efficient way of storing this data, will be to pre-compute the hourly fraction raster and store it in the database.

mwdunlap2004 commented 1 month ago

First Method, this one does the prism raster disaggregated to hourly and then resampled to our daymet resolution

This takes around 30 seconds all together, the ability to use our daily daymet data as well, means that we don't have the full amount stored, it's just the single daily raster.

We've done this step before, it just multiplies our prism daily raster by nldas2 to get rain for every hour

create temp table tmp_prism_fraction_raster as SELECT tstime, tsendtime,st_mapalgebra( ST_Resample( st_clip( fraction_table.rast, st_envelope(daily_table.reclass_rast), TRUE), daily_table.rast ),1, ST_SetBandNoDataValue(daily_table.reclass_rast,-9999), 1, '[rast1] * [rast2]' ) as rast from tmp_fraction_raster as fraction_table left outer join tmp_prism_feb12_daily as daily_table on (1=1);

This resamples our prism raster using our daymet daily raster as its reference, this seems to be quite a bit slower than when we have this step first, which makes sense as we have to do 24 times as much work.

create temp table tmp_prism_resample as select prism.tstime, prism.tsendtime, ST_Resample(prism.rast, reference.rast), prism.rast as rast from tmp_prism_fraction_raster as prism,tmp_daymet_feb12_daily as reference;

mwdunlap2004 commented 1 month ago

This method seems to be much faster, around 10-15 seconds for all the work compared to 30 I think not having to do 24 rasters at the end is just going to be more efficient, overall I think this is the method we should take.

This resamples our prism raster using our daily daymet raster as its reference, this runs extremely fast, finishing in just a few seconds. It creates one raster at the end, which we then disaggregate into hourly in the next step.

create temp table tmp_prism_resample as select prism.tstime, prism.tsendtime, ST_Resample(prism.rast, reference.rast), prism.reclass_rast as rast from tmp_prism_feb12_daily as prism, tmp_daymet_feb12_daily as reference;

create temp table tmp_prism_fraction_raster as SELECT tstime, tsendtime,st_mapalgebra( ST_Resample( st_clip( fraction_table.rast, st_envelope(daily_table.rast), TRUE), daily_table.rast ),1, ST_SetBandNoDataValue(daily_table.rast,-9999), 1, '[rast1] * [rast2]' ) as rast from tmp_fraction_raster as fraction_table left outer join tmp_prism_resample as daily_table on (1=1);

rburghol commented 1 month ago

@mwdunlap2004 that sounds good, let me understand, am I correct that:

mwdunlap2004 commented 1 month ago

@rburghol I think most of your statements are accurate, and I did have a typo in my second comment, in both examples we are using daymet as our reference raster. I will say, even though both are examples of prism, I do think this would apply to daymet, namely, it seems like it would be faster (if if was beneficial) to resample daymet and then do our multiplication to hourly as opposed to the other way around.

mwdunlap2004 commented 1 month ago

@rburghol These are the statistics for the four tables made in this process, the what I'm learning from this is that we cannot seem to sum the values and get what the original amount is (in this case I'm using our prism_daily_feb12, and the resampled prism daily as the base) When I sum the columns for the other tables we don't get the same value. My thought is the question we had before about what to do if an nldas2 column is 0 and we multiply it by PRISM with data. I think this is what is happening, some parts of nldas2 don't have rain maybe for an hour, or the whole day but PRISM does, which means we are losing data. The reason I believe this, is our resample of prism to daymet's resolution has the same values as our original prism raster, which means that the change in the data has to be happening in the multiplication process.

Screenshot 2024-08-07 at 8 34 41 AM Screenshot 2024-08-07 at 8 33 44 AM Screenshot 2024-08-07 at 8 34 58 AM

After we've resampled daily prism to daymet's resolution

These numbers are the same as our original prism_feb12_daily raster, which is a good sign.

Screenshot 2024-08-07 at 8 35 14 AM
mwdunlap2004 commented 1 month ago

This is an attempt at running these steps on a single watershed, for the purposes of showing our progress and creating a couple of rasters, The set up steps are located in #71

We've done this step before, it just multiplies our prism daily raster by nldas2 to get rain for every hour

create temp table tmp_prism_fraction_01665500_raster as
SELECT tid, tstime, tsendtime,st_mapalgebra(
  ST_Resample(
    st_clip(
      fraction_table.rast,
      st_envelope(daily_table.rast),
    TRUE), 
    daily_table.rast
  ),1, 
  ST_SetBandNoDataValue(daily_table.rast,-9999), 1, 
  '[rast1] * [rast2]'
) as rast
from tmp_nldas2_fraction_raster_01665500 as fraction_table
left outer join tmp_prism_01665500 as daily_table
on (1=1);

This resamples our prism raster using our daymet daily raster as its reference, this seems to be quite a bit slower than when we have this step first, which makes sense as we have to do 24 times as much work.

create temp table tmp_prism_resample_01665500 as
select prism.tstime, prism.tsendtime, ST_Resample(prism.rast, reference.rast), prism.rast as rast
from tmp_prism_fraction_01665500_raster as prism,tmp_daymet_01665500 as reference;
create temp table tmp_daymet_fraction_01665500_raster as
SELECT tid, tstime, tsendtime,st_mapalgebra(
  ST_Resample(
    st_clip(
      fraction_table.rast,
      st_envelope(daily_table.rast),
    TRUE), 
    daily_table.rast
  ),1, 
  ST_SetBandNoDataValue(daily_table.rast,-9999), 1, 
  '[rast1] * [rast2]'
) as rast
from tmp_nldas2_fraction_raster_01665500 as fraction_table
left outer join tmp_daymet_01665500 as daily_table
on (1=1);