Gauge adjustment

[parichat]

Hi all,

I used adjustment module of wradlib to correct bias of rain radar field.
The radar rain is average from 6-minute time interval of radar operation into hourly (mm/hr) and the rain gauges also in hourly.
I think the rain gauge network is quite dense enough, 442 rain gauges within radar coverage area 480 km x 480 km.

First, I tried on "Mean field bias (MFB)" but it shows not much improvement (attached file: MFB.png).

Second tried on "Additive method", it shows perfect result. Then I validated the result field by independent rain gauge, it show worse correlation than before adjustment (attached file: add.png).

Do you have any idea?  Suggest me please?

Best regards,
Parichat

[Maik Heistermann]

Dear Parichat,

The additive method has the property to perfectly reproduce the rain gage observations at the rain gage locations. If you compare against the rain gage data you
used for adjustment, you will get a perfect agreement. So you did the right thing and compared the adjusted field against independent rain gage observations. However, it is
hard to evaluate the quality of the adjustment from just one time step. The MFB adjustment shows that your radar not only suffers from severe underestimation (bias), but that the error also has a strong random component (after bias correction by MFB, you see a massive scatter). In such a situation, every adjustment method will be challenged. On the other hand, it also might be the case that an additive error model for radar correction is not adequate in this particular case. Maybe a multiplicative error model (variable in space) will be more adequate. In general, an additive error model is not very well behaved as it changes the distribution function because negative values from adjustment have to be truncated.

I have to admit that the AdjustAdd class is not yet thoroughly tested. The adjust module is still experimental...so there might be bugs still. I will try to look into the AdjustAdd class as soon as possible in order to check for bugs. Beyond, I will try to implement a class AdjustMultiply so you can check the efffectiveness of a multiplicative error model. However - as said before - you should be careful to evaluate the success of an adjustment from only one time step.

Best regards,
Maik

[parichat]

Dear Maik,


Thank you for your useful reply.
I understood that one time step couldn't say it successful.

So I did more long term analysis by expanding resolution of radar from 1x1 km to 4x4 km grid.
And accumulate the rain radar from 6-minute time step to hourly, then compared again for both methods.
And the result of scatter plot as in picture.

How do you think about it?

P.S. I'm waiting "AdjustMultiply" module and after launch I will test it for you.


Best regards,
Parichat

[Maik Heistermann]

Dear Parichat,

thanks for providing your new results! I am non 100% sure if I
understand your figure correctly: I guess the four rows represent four
consequtive days, from Sep 22-25 2009? Does every scatter plot
represent 24 hourly time steps, thus (24 x number of gages)?

Unfortunately, the resolution of the image is not high enough to read
the quality measures (correlation, RMSE, ...). From a pure visual
inspection, I would say, that the results are not too good, but that I
have seen worse. At least, the adjusted results seem to be a bit
better than the original ones. However, four days still are quite
short, typical verification studies span at least several years of
data ;-)

Some additonal remarks:

1. Why did you resample your data to 4 km2 resolution?

2. There seems to be a systematic underestimation of high
intensities...this can have many reasons, one could be an inadequate
Z-R- relationship, another can be attenuation by heavy rain. Do you
use C-band or S-band radar? For C-band you have to expect that
attenuation plays a major role in case of high intensities. In that
case, you could try the wradlib.atten module for attenuation correction.

3. You can try to play around with the parameters of the adjustment
approaches, e.g. the invserse distance power for the interpolation of
the error etc. (see library reference). Maybe you can imporve your
results.

4. I assume that the scatter plots under the "Orginal" column are
created using the gages that you used for adjustment while the
"Validation" column contains scatter plots created from the validation
gages. This makes it hard to evaluate the success of the adjustment.
In order to evaluate the success, the quality measures for "original"
and "validation" should be computed from using only one - the
validation - set of gages.

Finally, I pushed quite some changesets yesterday, including some
changes in the additive method (in case invalid values are present)
and, more important, the multiplicative adjustment via AdjustMultiply.
So please get the latest version of wradlib from
http://bitbucket.org/wradlib/wradlib/get/6faf016177c9.zip. Maybe you
can also try AdjusMultiply and report your results?

Thanks a lot, again, for sharing your experiences with the adjustment module.

Best regards,
Maik

> --
> You received this message because you are subscribed to the Google
> Groups "wradlib-users" group.
> To unsubscribe from this group, send email to
> wradlib-user...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>



--
Dr. Maik Heistermann
Universität Potsdam
Institut für Erd- und Umweltwissenschaften
Karl-Liebknecht-Str. 24-25
14476 Potsdam-Golm
phone: +49 331 977 2671

[parichat]

Dear Maik,

>>  1. Why did you resample your data to 4 km2 resolution?

          Because I would like to reduce the spatial distribution error.

>> 2. There seems to be a systematic underestimation of high  
intensities...this can have many reasons, one could be an inadequate  
Z-R- relationship, another can be attenuation by heavy rain. Do you  
use C-band or S-band radar? For C-band you have to expect that  
attenuation plays a major role in case of high intensities. In that  
case, you could try the wradlib.atten module for attenuation correction.

         Yes, you are right. It is underestimation at high intensity, I have no idea about that.
For my radar data is S-band, I try to use wradlib.atten but results show some large value such as "1.798e+308" as in the attached file "atten.txt".
My input data are in first column and results in final column, the used this command "k = atten.correctAttenuationConstrained(data, mode='zero')".

Did I do something wrong?

>> 3. You can try to play around with the parameters of the adjustment  
approaches, e.g. the invserse distance power for the interpolation of  
the error etc. (see library reference). Maybe you can imporve your  
results.

    I have been tried on one case for sensitivity test of p value as show in the attached file, "p.png", and found that best p is 2 which give less Root Mean Square Error (RMSE).

>> 4. I assume that the scatter plots under the "Orginal" column are  
created using the gages that you used for adjustment while the  
"Validation" column contains scatter plots created from the validation  
gages. This makes it hard to evaluate the success of the adjustment.  
In order to evaluate the success, the quality measures for "original"  
and "validation" should be computed from using only one - the  
validation - set of gages.

      Yes, your understood is correct. I applied those three method for a extreme rainfall continuous 3 days and found Multiplication gave best results among them.
I pick up an hourly to show you as in picture name add, mfb and multi in attached file.
      How do you think?
      And I have a question, why the original (means radar match up gauges) have good correlation than others even through didn't adjust yet??


Best regards,
Parichat

[Maik Heistermann]

Dear Parichat,

thanks for the new results - I think they're not too bad!

Just some remarks:

1. I think the resampling to a coarser resolution is not necessary
because the adjustment methods computes the radar value at the gage
locations based on a set of neighbouring radar bins. You can also
increase the number of neighbours if you want to achieve a smoothing
effect in space (default value for argument nnear_raws: 9) - please
see library reference.

2. Maybe you could try another Z-R relation, e.g. Z = 300 R^1.4 as
used by NOOA to reduce ovserestimation of light rainfall and
underestimation of heavy rainfall, or NOOA's Z-R relation for tropical
rainfall conditions (Z = 250 R^1.2). However, many of the "simple" Z-R
related effects should be accounted for in the course of adjustment.

3. I don't know about what went wrong in the attenuation correction,
however, since you are using an S-Band radar, attenuation by heavy
rainfall can be considered as negligible. So I would drop the
attenuation correction anyways. Could you maybe tell me the location
and the operator of your radar - just out of curiosity?

4. I'd like to stress again that in order to evaluate the success of
your adjustment, you should compute the performance measures for
unadjusted and adjusted radar based on the same set of validation rain
gages. You can also use the wradlib.verify module for computing the
error metrics (wradlib.verify.ErrorMetrics), but you also selected
useful metrics for your plots, anyway. The inconsistent correlation
could be a result of the different rain gage sets for adjustment and
validation.

5. I think the example you provided for AdjustMultiply looks quite
promising. Of course, it is hard to correct those radar values which
are almost zero by a multiplicative approach. Which value did you use
for the function argument minval in AdjustMultiply.__call__()? This
threshold is intended to avoid unrealistic ratio values caused by very
low radar values. For such situations as yours, a combination of
additive and multiplicative adjustment would be helpful. I am working
on implementing a prototype, but I can't promise a deadline, yet.

Best regards,
Maik



Quoting parichat <pari...@caos-a.geophys.tohoku.ac.jp>:

[Thomas Pfaff]

Hi Parichat,

 

I also don’t know, where those large values come from (might be a bug of some sort), but I would like to make two comments for attenuation correction in general.

 

1)       for S-Band, attenuation has been shown to be negligible in most rainfall situations, so you might not make a big mistake by ignoring it altogether

2)       all functions in the atten module implicitly assume you are using C-Band data, because the default parameters for the k-Z relation are those we (or better Stefan Krämer) found for C-Band.
As the parameters are highly dependent on the wavelength, they will produce nonsense results for your S-Band data. As long as you don’t find parameters for S-Band, you should not use these functions.

 

Have you checked the distance between the radar and the gauges?

Depending on the maximum observation range of the radar, the elevation angle and the distance to the gauges, the underestimates may just as well be due to vertical profile effects.

 

I also had a look at your sensitivity analysis for the interpolation power parameter p.

I would be interested how the results would look like for powers smaller than 2. Maybe they would get even better?

 

Best regards

 

 

Thomas

 

---

[Maik Heistermann]

Dear Parichat,

I just pushed a new adjustment method called AdjustMixed to the wradlib repo (see http://wradlib.bitbucket.org/adjust.html). Feel free to test it. AdjustMixed assumes that you have both additive and multiplicative error terms and tries to optimise their combination. Again, it would be great if you could post the results in case you use that approach. Good luck!

Best regards,
Maik

[parichat]

Dear Maik and Thomas,

So sorry for late reply.
I don't know why my post on 2 day ago was deleted twice, then I give up the post.

1.   I used a coarser resolution because I will combine the rainfall from radar to satellite rainfall which resolution is 4x4 km grid.

2.&4.  The attached picture shows verification by wradlib.verify.ErrorMetrics as your require.
Actually, I classified cloud into 2 types, convective and stratiform by a method which provide by Steiner (1995) as result in file name "classification.jpg".
Then use Z-R relationship to convert dBz to rainfall by separately parameter specific for each type as in picture "verify.jpg".

3.  My radar located at Phimai, Thailand (15º10’90’’N, 102º33’88’’E).

5. For AdjustMultiply threshol value, I use 0.5 because it is threshold of rain gauge.


I will try to examine the correlation of distance between gauge  and radar, also the sensitivity of p power as Thomas recommends.
And I promise that I will show result of new adjust method which you just added.

 
Thanks.
Parichat

[Maik Heistermann]

Dear Parichat,

thanks a lot for your latest results. It is really great that you keep
us so well informed about your experiences and your progress - thanks!

The verification results look VERY promising. Are these actually the
results for the independent verification gages or are these based on
the gages that you used for adjustment?

Just out of curiosity (only answer if you want to and if you have the time!):

1. Which satellite product do you intend to use?

2. What is your general strategy for combining radar and satellite products?

And one more thing: Would you be willing to share your code fragments
implementing the Steiner method for stratiform/convective separation?
If yes, we could implement this functionality into wradlib - which
would be really useful.

Good luck!

Maik

> --
> You received this message because you are subscribed to the Google
> Groups "wradlib-users" group.
> To unsubscribe from this group, send email to
> wradlib-user...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>
>

[parichat]

1. & 2. I will use Brightness temperature from MTSAT geostationary satellite and convert them into rain rate similar as radar.  I could not say deep detail because it is not success yes, still ongoing.

I wanna share my code but it wrote in fortran because I just started to use python after I found wradlib. 


Best,
Parichat

[Maik Heistermann]

Dear Parichat,

no problem, we can translate the code to Python.

Bests,
Maik

>> > wradlib-user...@googlegroups.com <javascript:>.

[parichat]

Dear Maik,

The attached file is my cloud classification code.
By the way, "AdjustMixed" module couldn't work normally.
It shows error message as following:

=============================
Traceback (most recent call last):
  File "adjust.py", line 101, in <module>
    adjuster = adjust.AdjustMixed(g_coords, r_coords, r_coords, stat='best', p_idw=2., mingages=5)
  File "/Library/Python/2.6/site-packages/wradlib/adjust.py", line 82, in __init__
    self.get_raw_at_obs = Raw_at_obs(self.obs_coords,  self.raw_coords, nnear=nnear_raws, stat=stat)
  File "/Library/Python/2.6/site-packages/wradlib/adjust.py", line 606, in __init__
    self.raw_ix = _get_neighbours_ix(obs_coords, raw_coords, nnear)
  File "/Library/Python/2.6/site-packages/wradlib/adjust.py", line 645, in _get_neighbours_ix
    return tree.query(obs_coords, k=nnear)[1]
  File "ckdtree.pyx", line 514, in scipy.spatial.ckdtree.cKDTree.query (scipy/spatial/ckdtree.c:4529)
TypeError: only length-1 arrays can be converted to Python scalars
===============================

Best,
Parichat

[Maik Heistermann]

Dear Parichat,

thanks a lot for the code.

Concerning the error you found: I am a bit surprised that you get an
error on the initialisation of the AdjustMixed instance as it is
initialised in the very same way as all the other adjustment classes.
I guess the error is in the following line you used:

adjuster = adjust.AdjustMixed(g_coords, r_coords, r_coords,
stat='best', p_idw=2., mingages=5)

In this line, you pass r_coords - the array of radar acoordinates -
twice. Instead of the second r_coords, you should pass the argument
nnear_raw which is the number of nearest neighbours in the radar grid
which is used to compute the radar values at the gage location
(deafult value: 9).

I hope this helps.

[parichat]

Oh! I'm sorry that I didn't check carefully.

By the way, if I would like to plot distance of gauge from radar on x axis and correlation on y axis.
You library have any module can do it?

And I saw the "ver.py" to plot vertical profile, is it available for now?
Because I wanna plot cross section to see the data in vertical.

Best,
Parichat