Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Reversibility of Spatial Resample and Temporal Disaggregate Steps #79

Open
2 tasks
rburghol opened this issue Aug 5, 2024 · 6 comments
Open
2 tasks
Assignees

Comments

@rburghol
Copy link
Contributor

rburghol commented Aug 5, 2024

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.

  • Dissag/Resample workflow tests:
    • compare Dissag PRISM to hourly -> Resample Disagg PRISM to 1km x 1km" to "Resample PRISM to 1km x 1km -> Dissag Resampled PRISM to Hourly
@mwdunlap2004
Copy link

mwdunlap2004 commented Aug 6, 2024

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
Copy link

mwdunlap2004 commented Aug 6, 2024

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
Copy link
Contributor Author

rburghol commented Aug 6, 2024

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

  • For PRISM, we are (and rightly so) re-sampling the NLDAS2 fraction table to the resolution of the PRISM daily table.
  • This PRISM resample is faster than doing it for daymet, which also makes sense since daymet has 16 times the number of cells per unit area as PRISM.
  • We will still have to do this process with daymet, so that will be slow, but that's just the price we have to pay :),
  • Question still remains, however, does our disaggregation and resembling process order result in different values?

@mwdunlap2004
Copy link

mwdunlap2004 commented Aug 7, 2024

@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.

  • For both examples we are resampling the nldas2 fraction table to the resolution of the prism daily table, in the first example this is at prisms original resolution (4km x 4km), in the second example we've resampled to daymet's resolution (I had a typo in the explanation but not the code itself).
  • The resample to daymet first is faster, since we don't have to resample 24 rasters to that resolution we are only doing 1, and then doing multiplication and clipping which I guess is faster. The resample takes a couple seconds and the multiplication takes around 10 in the second example.
  • I'm not sure if we have to do this with daymet, I know we must do the multiplication and resample of nldas2, but if we are using daymet's resolution as our reference then I would assume that means we didn't have to.
  • I've done a little bit of statical analysis, see my next comment, once I am able to output our rasters so we can see them I'll be able to do more.

@mwdunlap2004
Copy link

@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
Copy link

mwdunlap2004 commented Aug 8, 2024

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);

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants