[Met_help] [rt.rap.ucar.edu #46080] History for Vertical Interpolation and MADIS Observations

RAL HelpDesk {for John Halley Gotway} met_help at ucar.edu
Tue Apr 19 08:52:42 MDT 2011


----------------------------------------------------------------
  Initial Request
----------------------------------------------------------------

Hi,

As you can probably tell by the subject, I have a few questions to ask 
about MET.  I'll start off with my shorter question...

I was wondering whether or not anyone at DTC has investigated using the 
MADIS data stream as a point observation source.  The data come through 
our feed already in netCDF, but they're not in the correct netCDF format 
to use them directly.  I figured the best approach would be to write 
them out in ascii and then use the ASCII2NC tool to convert them back to 
netCDF.  I could write something to do this myself, but I figured it 
wouldn't hurt to ask if anyone out there has any NCL scripts or 
something of that nature to convert these data.

My second question involves how the vertical interpolation is done in 
point_stat.  In the old version of MET (very old, i.e. v0.9), it was 
easy to see in the code that if you have an observation with a pressure 
between two pressure levels, the forecast fields would be interpolated 
to the level of the observation.  However, if an observation wasn't 
bounded by two pressure levels, it would just be compared directly to 
the nearest pressure level.  For example, if your grib file contains 
1000 hPa as the lowest pressure level, and an observation has a pressure 
of 1020 hPa, those two would be directly compared without any 
interpolation/extrapolation.

This also seems to be the case in METv3.0.  I see it says...

/// Check the levels for the fcst and obs fields.  If the
       // fcst_field is a range of pressure levels, check to see if the
       // range of obs_field pressure levels is wholly contained in the
       // fcst levels.  If not, print a warning message.
       if(gc_pd[i].fcst_gci.lvl_type == PresLevel &&
          gc_pd[i].fcst_gci.lvl_1    != gc_pd[i].fcst_gci.lvl_2 &&
          (gc_pd[i].obs_gci.lvl_1 <  gc_pd[i].fcst_gci.lvl_1 ||
           gc_pd[i].obs_gci.lvl_2 >  gc_pd[i].fcst_gci.lvl_2)) {

          cout << "***WARNING***: PointStatConfInfo::process_config() -> "
<< "The range of requested observation pressure levels "
<< "is not contained within the range of requested "
<< "forecast pressure levels.  No vertical interpolation "
<< "will be performed for observations falling outside "
<< "the range of forecast levels.  Instead, they will be "
<< "matched to the single nearest forecast level.\n\n"
<< flush;

/I don't see exactly where the calculation is being done in the code to 
prove or disprove this.  When I ran a comparison between the two 
versions of MET using the same obs and forecast files, I used ADPSFC 
observations and TMP/P1000-700, which is what you needed to do in the 
earlier version to make sure you weren't getting rid of surface 
observations at elevation.  With METv0.9, there were 516 observations 
used within my domain of interest, compared to 181 observations with 
METv3.0.  So it seems like something else is going on that isn't 
immediately clear.  Ideally, I would just use the TMP/Z2, 2-m diagnosed 
forecast field to compare to the ADPSFC observations (this method yields 
690 observations) instead of the P1000-700 method, but not all the model 
products I use have that diagnosed field.  Anyway....if you could 
clarify what is actually being done with the vertical interpolation, it 
would definitely put my mind at ease.

Thank you again for all your help.


- Jeff Zielonka


----------------------------------------------------------------
  Complete Ticket History
----------------------------------------------------------------

Subject: Re: [rt.rap.ucar.edu #46080] Vertical Interpolation and MADIS Observations
From: John Halley Gotway
Time: Mon Apr 18 12:45:32 2011

Jeff,

Currently, no, we are not using MADIS observations within the DTC as a
point observation source.  However, I do know of a couple of other MET
users using the NetCDF MADIS observations.  In fact, one
user wrote a C++ program (madis2nc) to read in the MADIS observations
and write some of them out in a NetCDF format that MET can read.
Actually, another user requested that code last week.  So I got
some sample data from that user and am trying to run it through the
conversion.  But I've run into some problems with how the timing
information is stored in the input NetCDF file.

I'll bring this issue up to the MET developers group and let them know
that several users are requesting support for MADIS observations.  In
the shorter term, I'll let you know if we get this madis2nc
code working well and will be happy to send it to you for your own
testing.

Alternatively, writing an NCL script to process this data should work
too.

Regarding how matching in the vertical is done...

Most of the logic for how point observations are processed is
contained in METv3.0.1/lib/vx_met_util/pair_data.cc in the routine:

   void GCPairData::add_obs(float *hdr_arr,     char *hdr_typ_str,
                            char  *hdr_sid_str, unixtime hdr_ut,
                            float *obs_arr,     Grid &gr)

In there, you'll see the following code for handling the level
information - including pressure levels, accumulation intervals,
vertical levels (e.g. meters AGL), and miscellaneous level types:

   // For pressure levels, check if the observation pressure level
   // falls in the requsted range.
   if(obs_gci.lvl_type == PresLevel) {

      if(obs_lvl < obs_gci.lvl_1 ||
         obs_lvl > obs_gci.lvl_2) {
         rej_lvl++;
         return;
      }
   }
   // For accumulations, check if the observation accumulation
interval
   // matches the requested interval.
   else if(obs_gci.lvl_type == AccumLevel) {

      if(obs_lvl < obs_gci.lvl_1 ||
         obs_lvl > obs_gci.lvl_2) {
         rej_lvl++;
         return;
      }
   }
   // For vertical levels, check for a surface message type or if the
   // observation height falls within the requested range.
   else if(obs_gci.lvl_type == VertLevel) {

      if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL &&
         (obs_hgt < obs_gci.lvl_1 ||
          obs_hgt > obs_gci.lvl_2)) {
         rej_lvl++;
         return;
      }
   }
   // For all other level types (RecNumber, NoLevel), check for a
   // surface message type or if the observation height falls within
   // the requested range.
   else {

      if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL &&
         (obs_hgt < obs_gci.lvl_1 ||
          obs_hgt > obs_gci.lvl_2)) {
         rej_lvl++;
         return;
      }
   }

Hopefully that helps clarify what's currently being done in METv3.0.1.
As for the differences between METv0.9 and METv3.0.1, I really can't
say for sure.

Based on the logic listed there, I do think that if you have the
following settings...
  fcst_field[] = [ "TMP/P1000-700" ];
  obs_field[]  = [ "TMP/P1100-600" ];
And the code comes across an observation of temperature at 1050, it'd
just match it to the forecast value at 1000.

Yes, I do agree that it would be preferable to just verify 2-meter
temperature (TMP/Z2) against ADPSFC observations, but I understand
that you're limited by what model fields you have available.

If it would be helpful to you, feel free to pick a case where you get
different numbers of matched pairs from METv0.9 and METv3.0.1.  Send
me a single forecast GRIB file, the corresponding point
observation file, and config files for METv0.9 and METv3.0.1.  I could
run it here and figure out why the numbers differ.  If you'd like to
do that, you can post data to our ftp site using the
following instructions:
  http://www.dtcenter.org/met/users/support/met_help.php#ftp

Thanks,
John Halley Gotway
met_help at ucar.edu


On 04/13/2011 07:56 AM, RAL HelpDesk {for Jeff Zielonka} wrote:
>
> Wed Apr 13 07:56:25 2011: Request 46080 was acted upon.
> Transaction: Ticket created by zielonka at meteo.psu.edu
>        Queue: met_help
>      Subject: Vertical Interpolation and MADIS Observations
>        Owner: Nobody
>   Requestors: zielonka at meteo.psu.edu
>       Status: new
>  Ticket <URL:
https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080 >
>
>
> Hi,
>
> As you can probably tell by the subject, I have a few questions to
ask
> about MET.  I'll start off with my shorter question...
>
> I was wondering whether or not anyone at DTC has investigated using
the
> MADIS data stream as a point observation source.  The data come
through
> our feed already in netCDF, but they're not in the correct netCDF
format
> to use them directly.  I figured the best approach would be to write
> them out in ascii and then use the ASCII2NC tool to convert them
back to
> netCDF.  I could write something to do this myself, but I figured it
> wouldn't hurt to ask if anyone out there has any NCL scripts or
> something of that nature to convert these data.
>
> My second question involves how the vertical interpolation is done
in
> point_stat.  In the old version of MET (very old, i.e. v0.9), it was
> easy to see in the code that if you have an observation with a
pressure
> between two pressure levels, the forecast fields would be
interpolated
> to the level of the observation.  However, if an observation wasn't
> bounded by two pressure levels, it would just be compared directly
to
> the nearest pressure level.  For example, if your grib file contains
> 1000 hPa as the lowest pressure level, and an observation has a
pressure
> of 1020 hPa, those two would be directly compared without any
> interpolation/extrapolation.
>
> This also seems to be the case in METv3.0.  I see it says...
>
> /// Check the levels for the fcst and obs fields.  If the
>        // fcst_field is a range of pressure levels, check to see if
the
>        // range of obs_field pressure levels is wholly contained in
the
>        // fcst levels.  If not, print a warning message.
>        if(gc_pd[i].fcst_gci.lvl_type == PresLevel &&
>           gc_pd[i].fcst_gci.lvl_1    != gc_pd[i].fcst_gci.lvl_2 &&
>           (gc_pd[i].obs_gci.lvl_1 <  gc_pd[i].fcst_gci.lvl_1 ||
>            gc_pd[i].obs_gci.lvl_2 >  gc_pd[i].fcst_gci.lvl_2)) {
>
>           cout << "***WARNING***:
PointStatConfInfo::process_config() -> "
> << "The range of requested observation pressure levels "
> << "is not contained within the range of requested "
> << "forecast pressure levels.  No vertical interpolation "
> << "will be performed for observations falling outside "
> << "the range of forecast levels.  Instead, they will be "
> << "matched to the single nearest forecast level.\n\n"
> << flush;
>
> /I don't see exactly where the calculation is being done in the code
to
> prove or disprove this.  When I ran a comparison between the two
> versions of MET using the same obs and forecast files, I used ADPSFC
> observations and TMP/P1000-700, which is what you needed to do in
the
> earlier version to make sure you weren't getting rid of surface
> observations at elevation.  With METv0.9, there were 516
observations
> used within my domain of interest, compared to 181 observations with
> METv3.0.  So it seems like something else is going on that isn't
> immediately clear.  Ideally, I would just use the TMP/Z2, 2-m
diagnosed
> forecast field to compare to the ADPSFC observations (this method
yields
> 690 observations) instead of the P1000-700 method, but not all the
model
> products I use have that diagnosed field.  Anyway....if you could
> clarify what is actually being done with the vertical interpolation,
it
> would definitely put my mind at ease.
>
> Thank you again for all your help.
>
>
> - Jeff Zielonka

------------------------------------------------
Subject: Vertical Interpolation and MADIS Observations
From: John Halley Gotway
Time: Mon Apr 18 14:58:08 2011

Jeff,

We were able to resolve this issues with the madis2nc code from that
other user.  I was able to reformat some MADIS NetCDF data and use it
in a run of Point-Stat.

I looked closer at the madis2nc code and see that this really only
supports MADIS data that contains METARS.  Specifically, it looks for
this global attribute:
  title = "MADIS - Meteorological Surface - METAR"

If we add better support for this in a future release of MET, we'll
add in support for the other MADIS data types.

I've attached the a tar file containing the madis2nc code for
METv3.0.1.  Please perform the following steps:
  cd METv3.0.1
  mv /path/to/METv3.0.1_madis2nc_beta.tar.gz .
  tar -xvzf METv3.0.1_madis2nc_beta.tar.gz

And then try rebuilding MET.  The madis2nc executable should appear in
the bin directory.

Please note that I can't make any claims on the accuracy of this code.
It's contributed code at this point that we have spent little to no
time testing.  So please use it at your own risk.

Thanks and I hope this helps.

Just let me know if you have any questions.

Thanks,
John

On 04/18/2011 12:45 PM, RAL HelpDesk {for John Halley Gotway} wrote:
>
> <URL: https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080 >
>
> Jeff,
>
> Currently, no, we are not using MADIS observations within the DTC as
a point observation source.  However, I do know of a couple of other
MET users using the NetCDF MADIS observations.  In fact, one
> user wrote a C++ program (madis2nc) to read in the MADIS
observations and write some of them out in a NetCDF format that MET
can read.  Actually, another user requested that code last week.  So I
got
> some sample data from that user and am trying to run it through the
conversion.  But I've run into some problems with how the timing
information is stored in the input NetCDF file.
>
> I'll bring this issue up to the MET developers group and let them
know that several users are requesting support for MADIS observations.
In the shorter term, I'll let you know if we get this madis2nc
> code working well and will be happy to send it to you for your own
testing.
>
> Alternatively, writing an NCL script to process this data should
work too.
>
> Regarding how matching in the vertical is done...
>
> Most of the logic for how point observations are processed is
contained in METv3.0.1/lib/vx_met_util/pair_data.cc in the routine:
>
>    void GCPairData::add_obs(float *hdr_arr,     char *hdr_typ_str,
>                             char  *hdr_sid_str, unixtime hdr_ut,
>                             float *obs_arr,     Grid &gr)
>
> In there, you'll see the following code for handling the level
information - including pressure levels, accumulation intervals,
vertical levels (e.g. meters AGL), and miscellaneous level types:
>
>    // For pressure levels, check if the observation pressure level
>    // falls in the requsted range.
>    if(obs_gci.lvl_type == PresLevel) {
>
>       if(obs_lvl < obs_gci.lvl_1 ||
>          obs_lvl > obs_gci.lvl_2) {
>          rej_lvl++;
>          return;
>       }
>    }
>    // For accumulations, check if the observation accumulation
interval
>    // matches the requested interval.
>    else if(obs_gci.lvl_type == AccumLevel) {
>
>       if(obs_lvl < obs_gci.lvl_1 ||
>          obs_lvl > obs_gci.lvl_2) {
>          rej_lvl++;
>          return;
>       }
>    }
>    // For vertical levels, check for a surface message type or if
the
>    // observation height falls within the requested range.
>    else if(obs_gci.lvl_type == VertLevel) {
>
>       if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL &&
>          (obs_hgt < obs_gci.lvl_1 ||
>           obs_hgt > obs_gci.lvl_2)) {
>          rej_lvl++;
>          return;
>       }
>    }
>    // For all other level types (RecNumber, NoLevel), check for a
>    // surface message type or if the observation height falls within
>    // the requested range.
>    else {
>
>       if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL &&
>          (obs_hgt < obs_gci.lvl_1 ||
>           obs_hgt > obs_gci.lvl_2)) {
>          rej_lvl++;
>          return;
>       }
>    }
>
> Hopefully that helps clarify what's currently being done in
METv3.0.1.  As for the differences between METv0.9 and METv3.0.1, I
really can't say for sure.
>
> Based on the logic listed there, I do think that if you have the
following settings...
>   fcst_field[] = [ "TMP/P1000-700" ];
>   obs_field[]  = [ "TMP/P1100-600" ];
> And the code comes across an observation of temperature at 1050,
it'd just match it to the forecast value at 1000.
>
> Yes, I do agree that it would be preferable to just verify 2-meter
temperature (TMP/Z2) against ADPSFC observations, but I understand
that you're limited by what model fields you have available.
>
> If it would be helpful to you, feel free to pick a case where you
get different numbers of matched pairs from METv0.9 and METv3.0.1.
Send me a single forecast GRIB file, the corresponding point
> observation file, and config files for METv0.9 and METv3.0.1.  I
could run it here and figure out why the numbers differ.  If you'd
like to do that, you can post data to our ftp site using the
> following instructions:
>   http://www.dtcenter.org/met/users/support/met_help.php#ftp
>
> Thanks,
> John Halley Gotway
> met_help at ucar.edu
>
>
> On 04/13/2011 07:56 AM, RAL HelpDesk {for Jeff Zielonka} wrote:
>>
>> Wed Apr 13 07:56:25 2011: Request 46080 was acted upon.
>> Transaction: Ticket created by zielonka at meteo.psu.edu
>>        Queue: met_help
>>      Subject: Vertical Interpolation and MADIS Observations
>>        Owner: Nobody
>>   Requestors: zielonka at meteo.psu.edu
>>       Status: new
>>  Ticket <URL:
https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080 >
>>
>>
>> Hi,
>>
>> As you can probably tell by the subject, I have a few questions to
ask
>> about MET.  I'll start off with my shorter question...
>>
>> I was wondering whether or not anyone at DTC has investigated using
the
>> MADIS data stream as a point observation source.  The data come
through
>> our feed already in netCDF, but they're not in the correct netCDF
format
>> to use them directly.  I figured the best approach would be to
write
>> them out in ascii and then use the ASCII2NC tool to convert them
back to
>> netCDF.  I could write something to do this myself, but I figured
it
>> wouldn't hurt to ask if anyone out there has any NCL scripts or
>> something of that nature to convert these data.
>>
>> My second question involves how the vertical interpolation is done
in
>> point_stat.  In the old version of MET (very old, i.e. v0.9), it
was
>> easy to see in the code that if you have an observation with a
pressure
>> between two pressure levels, the forecast fields would be
interpolated
>> to the level of the observation.  However, if an observation wasn't
>> bounded by two pressure levels, it would just be compared directly
to
>> the nearest pressure level.  For example, if your grib file
contains
>> 1000 hPa as the lowest pressure level, and an observation has a
pressure
>> of 1020 hPa, those two would be directly compared without any
>> interpolation/extrapolation.
>>
>> This also seems to be the case in METv3.0.  I see it says...
>>
>> /// Check the levels for the fcst and obs fields.  If the
>>        // fcst_field is a range of pressure levels, check to see if
the
>>        // range of obs_field pressure levels is wholly contained in
the
>>        // fcst levels.  If not, print a warning message.
>>        if(gc_pd[i].fcst_gci.lvl_type == PresLevel &&
>>           gc_pd[i].fcst_gci.lvl_1    != gc_pd[i].fcst_gci.lvl_2 &&
>>           (gc_pd[i].obs_gci.lvl_1 <  gc_pd[i].fcst_gci.lvl_1 ||
>>            gc_pd[i].obs_gci.lvl_2 >  gc_pd[i].fcst_gci.lvl_2)) {
>>
>>           cout << "***WARNING***:
PointStatConfInfo::process_config() -> "
>> << "The range of requested observation pressure levels "
>> << "is not contained within the range of requested "
>> << "forecast pressure levels.  No vertical interpolation "
>> << "will be performed for observations falling outside "
>> << "the range of forecast levels.  Instead, they will be "
>> << "matched to the single nearest forecast level.\n\n"
>> << flush;
>>
>> /I don't see exactly where the calculation is being done in the
code to
>> prove or disprove this.  When I ran a comparison between the two
>> versions of MET using the same obs and forecast files, I used
ADPSFC
>> observations and TMP/P1000-700, which is what you needed to do in
the
>> earlier version to make sure you weren't getting rid of surface
>> observations at elevation.  With METv0.9, there were 516
observations
>> used within my domain of interest, compared to 181 observations
with
>> METv3.0.  So it seems like something else is going on that isn't
>> immediately clear.  Ideally, I would just use the TMP/Z2, 2-m
diagnosed
>> forecast field to compare to the ADPSFC observations (this method
yields
>> 690 observations) instead of the P1000-700 method, but not all the
model
>> products I use have that diagnosed field.  Anyway....if you could
>> clarify what is actually being done with the vertical
interpolation, it
>> would definitely put my mind at ease.
>>
>> Thank you again for all your help.
>>
>>
>> - Jeff Zielonka

------------------------------------------------
Subject: Re: [rt.rap.ucar.edu #46080] Vertical Interpolation and MADIS Observations
From: Jeff Zielonka
Time: Mon Apr 18 15:12:02 2011

Thanks John,

I'm tied up with a few things at the moment, but as soon as I get a
chance to put this in and test it, I'll let you know what I find.  I
did
just do a quick check, and "title : MADIS - Meteorological Surface -
METAR" does appear when I dump my obs files, so this has some promise.
I think the next most logical MADIS observation type to use would be
RAOB, and then more as users request them I suppose.  I'll also take a
look at the code and see if I can make any worthwhile updates.  Thanks
again.

- Jeff


On 4/18/2011 4:58 PM, RAL HelpDesk {for John Halley Gotway} wrote:
> Jeff,
>
> We were able to resolve this issues with the madis2nc code from that
other user.  I was able to reformat some MADIS NetCDF data and use it
in a run of Point-Stat.
>
> I looked closer at the madis2nc code and see that this really only
supports MADIS data that contains METARS.  Specifically, it looks for
this global attribute:
>    title = "MADIS - Meteorological Surface - METAR"
>
> If we add better support for this in a future release of MET, we'll
add in support for the other MADIS data types.
>
> I've attached the a tar file containing the madis2nc code for
METv3.0.1.  Please perform the following steps:
>    cd METv3.0.1
>    mv /path/to/METv3.0.1_madis2nc_beta.tar.gz .
>    tar -xvzf METv3.0.1_madis2nc_beta.tar.gz
>
> And then try rebuilding MET.  The madis2nc executable should appear
in the bin directory.
>
> Please note that I can't make any claims on the accuracy of this
code.  It's contributed code at this point that we have spent little
to no time testing.  So please use it at your own risk.
>
> Thanks and I hope this helps.
>
> Just let me know if you have any questions.
>
> Thanks,
> John
>
> On 04/18/2011 12:45 PM, RAL HelpDesk {for John Halley Gotway} wrote:
>> <URL: https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080>
>>
>> Jeff,
>>
>> Currently, no, we are not using MADIS observations within the DTC
as a point observation source.  However, I do know of a couple of
other MET users using the NetCDF MADIS observations.  In fact, one
>> user wrote a C++ program (madis2nc) to read in the MADIS
observations and write some of them out in a NetCDF format that MET
can read.  Actually, another user requested that code last week.  So I
got
>> some sample data from that user and am trying to run it through the
conversion.  But I've run into some problems with how the timing
information is stored in the input NetCDF file.
>>
>> I'll bring this issue up to the MET developers group and let them
know that several users are requesting support for MADIS observations.
In the shorter term, I'll let you know if we get this madis2nc
>> code working well and will be happy to send it to you for your own
testing.
>>
>> Alternatively, writing an NCL script to process this data should
work too.
>>
>> Regarding how matching in the vertical is done...
>>
>> Most of the logic for how point observations are processed is
contained in METv3.0.1/lib/vx_met_util/pair_data.cc in the routine:
>>
>>     void GCPairData::add_obs(float *hdr_arr,     char *hdr_typ_str,
>>                              char  *hdr_sid_str, unixtime hdr_ut,
>>                              float *obs_arr,     Grid&gr)
>>
>> In there, you'll see the following code for handling the level
information - including pressure levels, accumulation intervals,
vertical levels (e.g. meters AGL), and miscellaneous level types:
>>
>>     // For pressure levels, check if the observation pressure level
>>     // falls in the requsted range.
>>     if(obs_gci.lvl_type == PresLevel) {
>>
>>        if(obs_lvl<  obs_gci.lvl_1 ||
>>           obs_lvl>  obs_gci.lvl_2) {
>>           rej_lvl++;
>>           return;
>>        }
>>     }
>>     // For accumulations, check if the observation accumulation
interval
>>     // matches the requested interval.
>>     else if(obs_gci.lvl_type == AccumLevel) {
>>
>>        if(obs_lvl<  obs_gci.lvl_1 ||
>>           obs_lvl>  obs_gci.lvl_2) {
>>           rej_lvl++;
>>           return;
>>        }
>>     }
>>     // For vertical levels, check for a surface message type or if
the
>>     // observation height falls within the requested range.
>>     else if(obs_gci.lvl_type == VertLevel) {
>>
>>        if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL&&
>>           (obs_hgt<  obs_gci.lvl_1 ||
>>            obs_hgt>  obs_gci.lvl_2)) {
>>           rej_lvl++;
>>           return;
>>        }
>>     }
>>     // For all other level types (RecNumber, NoLevel), check for a
>>     // surface message type or if the observation height falls
within
>>     // the requested range.
>>     else {
>>
>>        if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL&&
>>           (obs_hgt<  obs_gci.lvl_1 ||
>>            obs_hgt>  obs_gci.lvl_2)) {
>>           rej_lvl++;
>>           return;
>>        }
>>     }
>>
>> Hopefully that helps clarify what's currently being done in
METv3.0.1.  As for the differences between METv0.9 and METv3.0.1, I
really can't say for sure.
>>
>> Based on the logic listed there, I do think that if you have the
following settings...
>>    fcst_field[] = [ "TMP/P1000-700" ];
>>    obs_field[]  = [ "TMP/P1100-600" ];
>> And the code comes across an observation of temperature at 1050,
it'd just match it to the forecast value at 1000.
>>
>> Yes, I do agree that it would be preferable to just verify 2-meter
temperature (TMP/Z2) against ADPSFC observations, but I understand
that you're limited by what model fields you have available.
>>
>> If it would be helpful to you, feel free to pick a case where you
get different numbers of matched pairs from METv0.9 and METv3.0.1.
Send me a single forecast GRIB file, the corresponding point
>> observation file, and config files for METv0.9 and METv3.0.1.  I
could run it here and figure out why the numbers differ.  If you'd
like to do that, you can post data to our ftp site using the
>> following instructions:
>>    http://www.dtcenter.org/met/users/support/met_help.php#ftp
>>
>> Thanks,
>> John Halley Gotway
>> met_help at ucar.edu
>>
>>
>> On 04/13/2011 07:56 AM, RAL HelpDesk {for Jeff Zielonka} wrote:
>>> Wed Apr 13 07:56:25 2011: Request 46080 was acted upon.
>>> Transaction: Ticket created by zielonka at meteo.psu.edu
>>>         Queue: met_help
>>>       Subject: Vertical Interpolation and MADIS Observations
>>>         Owner: Nobody
>>>    Requestors: zielonka at meteo.psu.edu
>>>        Status: new
>>>   Ticket<URL:
https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080>
>>>
>>>
>>> Hi,
>>>
>>> As you can probably tell by the subject, I have a few questions to
ask
>>> about MET.  I'll start off with my shorter question...
>>>
>>> I was wondering whether or not anyone at DTC has investigated
using the
>>> MADIS data stream as a point observation source.  The data come
through
>>> our feed already in netCDF, but they're not in the correct netCDF
format
>>> to use them directly.  I figured the best approach would be to
write
>>> them out in ascii and then use the ASCII2NC tool to convert them
back to
>>> netCDF.  I could write something to do this myself, but I figured
it
>>> wouldn't hurt to ask if anyone out there has any NCL scripts or
>>> something of that nature to convert these data.
>>>
>>> My second question involves how the vertical interpolation is done
in
>>> point_stat.  In the old version of MET (very old, i.e. v0.9), it
was
>>> easy to see in the code that if you have an observation with a
pressure
>>> between two pressure levels, the forecast fields would be
interpolated
>>> to the level of the observation.  However, if an observation
wasn't
>>> bounded by two pressure levels, it would just be compared directly
to
>>> the nearest pressure level.  For example, if your grib file
contains
>>> 1000 hPa as the lowest pressure level, and an observation has a
pressure
>>> of 1020 hPa, those two would be directly compared without any
>>> interpolation/extrapolation.
>>>
>>> This also seems to be the case in METv3.0.  I see it says...
>>>
>>> /// Check the levels for the fcst and obs fields.  If the
>>>         // fcst_field is a range of pressure levels, check to see
if the
>>>         // range of obs_field pressure levels is wholly contained
in the
>>>         // fcst levels.  If not, print a warning message.
>>>         if(gc_pd[i].fcst_gci.lvl_type == PresLevel&&
>>>            gc_pd[i].fcst_gci.lvl_1    != gc_pd[i].fcst_gci.lvl_2&&
>>>            (gc_pd[i].obs_gci.lvl_1<   gc_pd[i].fcst_gci.lvl_1 ||
>>>             gc_pd[i].obs_gci.lvl_2>   gc_pd[i].fcst_gci.lvl_2)) {
>>>
>>>            cout<<  "***WARNING***:
PointStatConfInfo::process_config() ->  "
>>> <<  "The range of requested observation pressure levels"
>>> <<  "is not contained within the range of requested"
>>> <<  "forecast pressure levels.  No vertical interpolation"
>>> <<  "will be performed for observations falling outside"
>>> <<  "the range of forecast levels.  Instead, they will be"
>>> <<  "matched to the single nearest forecast level.\n\n"
>>> <<  flush;
>>>
>>> /I don't see exactly where the calculation is being done in the
code to
>>> prove or disprove this.  When I ran a comparison between the two
>>> versions of MET using the same obs and forecast files, I used
ADPSFC
>>> observations and TMP/P1000-700, which is what you needed to do in
the
>>> earlier version to make sure you weren't getting rid of surface
>>> observations at elevation.  With METv0.9, there were 516
observations
>>> used within my domain of interest, compared to 181 observations
with
>>> METv3.0.  So it seems like something else is going on that isn't
>>> immediately clear.  Ideally, I would just use the TMP/Z2, 2-m
diagnosed
>>> forecast field to compare to the ADPSFC observations (this method
yields
>>> 690 observations) instead of the P1000-700 method, but not all the
model
>>> products I use have that diagnosed field.  Anyway....if you could
>>> clarify what is actually being done with the vertical
interpolation, it
>>> would definitely put my mind at ease.
>>>
>>> Thank you again for all your help.
>>>
>>>
>>> - Jeff Zielonka


------------------------------------------------
Subject: Re: [rt.rap.ucar.edu #46080] Vertical Interpolation and MADIS Observations
From: John Halley Gotway
Time: Tue Apr 19 08:52:27 2011

Jeff,

Sounds good.  I'll go ahead and resolve this MET-Help ticket for now.
But if more questions arise, just let us know.

Thanks,
John

On 04/18/2011 03:12 PM, RAL HelpDesk {for Jeff Zielonka} wrote:
>
> <URL: https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080 >
>
> Thanks John,
>
> I'm tied up with a few things at the moment, but as soon as I get a
> chance to put this in and test it, I'll let you know what I find.  I
did
> just do a quick check, and "title : MADIS - Meteorological Surface -
> METAR" does appear when I dump my obs files, so this has some
promise.
> I think the next most logical MADIS observation type to use would be
> RAOB, and then more as users request them I suppose.  I'll also take
a
> look at the code and see if I can make any worthwhile updates.
Thanks
> again.
>
> - Jeff
>
>
> On 4/18/2011 4:58 PM, RAL HelpDesk {for John Halley Gotway} wrote:
>> Jeff,
>>
>> We were able to resolve this issues with the madis2nc code from
that other user.  I was able to reformat some MADIS NetCDF data and
use it in a run of Point-Stat.
>>
>> I looked closer at the madis2nc code and see that this really only
supports MADIS data that contains METARS.  Specifically, it looks for
this global attribute:
>>    title = "MADIS - Meteorological Surface - METAR"
>>
>> If we add better support for this in a future release of MET, we'll
add in support for the other MADIS data types.
>>
>> I've attached the a tar file containing the madis2nc code for
METv3.0.1.  Please perform the following steps:
>>    cd METv3.0.1
>>    mv /path/to/METv3.0.1_madis2nc_beta.tar.gz .
>>    tar -xvzf METv3.0.1_madis2nc_beta.tar.gz
>>
>> And then try rebuilding MET.  The madis2nc executable should appear
in the bin directory.
>>
>> Please note that I can't make any claims on the accuracy of this
code.  It's contributed code at this point that we have spent little
to no time testing.  So please use it at your own risk.
>>
>> Thanks and I hope this helps.
>>
>> Just let me know if you have any questions.
>>
>> Thanks,
>> John
>>
>> On 04/18/2011 12:45 PM, RAL HelpDesk {for John Halley Gotway}
wrote:
>>> <URL: https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080>
>>>
>>> Jeff,
>>>
>>> Currently, no, we are not using MADIS observations within the DTC
as a point observation source.  However, I do know of a couple of
other MET users using the NetCDF MADIS observations.  In fact, one
>>> user wrote a C++ program (madis2nc) to read in the MADIS
observations and write some of them out in a NetCDF format that MET
can read.  Actually, another user requested that code last week.  So I
got
>>> some sample data from that user and am trying to run it through
the conversion.  But I've run into some problems with how the timing
information is stored in the input NetCDF file.
>>>
>>> I'll bring this issue up to the MET developers group and let them
know that several users are requesting support for MADIS observations.
In the shorter term, I'll let you know if we get this madis2nc
>>> code working well and will be happy to send it to you for your own
testing.
>>>
>>> Alternatively, writing an NCL script to process this data should
work too.
>>>
>>> Regarding how matching in the vertical is done...
>>>
>>> Most of the logic for how point observations are processed is
contained in METv3.0.1/lib/vx_met_util/pair_data.cc in the routine:
>>>
>>>     void GCPairData::add_obs(float *hdr_arr,     char
*hdr_typ_str,
>>>                              char  *hdr_sid_str, unixtime hdr_ut,
>>>                              float *obs_arr,     Grid&gr)
>>>
>>> In there, you'll see the following code for handling the level
information - including pressure levels, accumulation intervals,
vertical levels (e.g. meters AGL), and miscellaneous level types:
>>>
>>>     // For pressure levels, check if the observation pressure
level
>>>     // falls in the requsted range.
>>>     if(obs_gci.lvl_type == PresLevel) {
>>>
>>>        if(obs_lvl<  obs_gci.lvl_1 ||
>>>           obs_lvl>  obs_gci.lvl_2) {
>>>           rej_lvl++;
>>>           return;
>>>        }
>>>     }
>>>     // For accumulations, check if the observation accumulation
interval
>>>     // matches the requested interval.
>>>     else if(obs_gci.lvl_type == AccumLevel) {
>>>
>>>        if(obs_lvl<  obs_gci.lvl_1 ||
>>>           obs_lvl>  obs_gci.lvl_2) {
>>>           rej_lvl++;
>>>           return;
>>>        }
>>>     }
>>>     // For vertical levels, check for a surface message type or if
the
>>>     // observation height falls within the requested range.
>>>     else if(obs_gci.lvl_type == VertLevel) {
>>>
>>>        if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL&&
>>>           (obs_hgt<  obs_gci.lvl_1 ||
>>>            obs_hgt>  obs_gci.lvl_2)) {
>>>           rej_lvl++;
>>>           return;
>>>        }
>>>     }
>>>     // For all other level types (RecNumber, NoLevel), check for a
>>>     // surface message type or if the observation height falls
within
>>>     // the requested range.
>>>     else {
>>>
>>>        if(strstr(onlysf_msg_typ_str, hdr_typ_str) == NULL&&
>>>           (obs_hgt<  obs_gci.lvl_1 ||
>>>            obs_hgt>  obs_gci.lvl_2)) {
>>>           rej_lvl++;
>>>           return;
>>>        }
>>>     }
>>>
>>> Hopefully that helps clarify what's currently being done in
METv3.0.1.  As for the differences between METv0.9 and METv3.0.1, I
really can't say for sure.
>>>
>>> Based on the logic listed there, I do think that if you have the
following settings...
>>>    fcst_field[] = [ "TMP/P1000-700" ];
>>>    obs_field[]  = [ "TMP/P1100-600" ];
>>> And the code comes across an observation of temperature at 1050,
it'd just match it to the forecast value at 1000.
>>>
>>> Yes, I do agree that it would be preferable to just verify 2-meter
temperature (TMP/Z2) against ADPSFC observations, but I understand
that you're limited by what model fields you have available.
>>>
>>> If it would be helpful to you, feel free to pick a case where you
get different numbers of matched pairs from METv0.9 and METv3.0.1.
Send me a single forecast GRIB file, the corresponding point
>>> observation file, and config files for METv0.9 and METv3.0.1.  I
could run it here and figure out why the numbers differ.  If you'd
like to do that, you can post data to our ftp site using the
>>> following instructions:
>>>    http://www.dtcenter.org/met/users/support/met_help.php#ftp
>>>
>>> Thanks,
>>> John Halley Gotway
>>> met_help at ucar.edu
>>>
>>>
>>> On 04/13/2011 07:56 AM, RAL HelpDesk {for Jeff Zielonka} wrote:
>>>> Wed Apr 13 07:56:25 2011: Request 46080 was acted upon.
>>>> Transaction: Ticket created by zielonka at meteo.psu.edu
>>>>         Queue: met_help
>>>>       Subject: Vertical Interpolation and MADIS Observations
>>>>         Owner: Nobody
>>>>    Requestors: zielonka at meteo.psu.edu
>>>>        Status: new
>>>>   Ticket<URL:
https://rt.rap.ucar.edu/rt/Ticket/Display.html?id=46080>
>>>>
>>>>
>>>> Hi,
>>>>
>>>> As you can probably tell by the subject, I have a few questions
to ask
>>>> about MET.  I'll start off with my shorter question...
>>>>
>>>> I was wondering whether or not anyone at DTC has investigated
using the
>>>> MADIS data stream as a point observation source.  The data come
through
>>>> our feed already in netCDF, but they're not in the correct netCDF
format
>>>> to use them directly.  I figured the best approach would be to
write
>>>> them out in ascii and then use the ASCII2NC tool to convert them
back to
>>>> netCDF.  I could write something to do this myself, but I figured
it
>>>> wouldn't hurt to ask if anyone out there has any NCL scripts or
>>>> something of that nature to convert these data.
>>>>
>>>> My second question involves how the vertical interpolation is
done in
>>>> point_stat.  In the old version of MET (very old, i.e. v0.9), it
was
>>>> easy to see in the code that if you have an observation with a
pressure
>>>> between two pressure levels, the forecast fields would be
interpolated
>>>> to the level of the observation.  However, if an observation
wasn't
>>>> bounded by two pressure levels, it would just be compared
directly to
>>>> the nearest pressure level.  For example, if your grib file
contains
>>>> 1000 hPa as the lowest pressure level, and an observation has a
pressure
>>>> of 1020 hPa, those two would be directly compared without any
>>>> interpolation/extrapolation.
>>>>
>>>> This also seems to be the case in METv3.0.  I see it says...
>>>>
>>>> /// Check the levels for the fcst and obs fields.  If the
>>>>         // fcst_field is a range of pressure levels, check to see
if the
>>>>         // range of obs_field pressure levels is wholly contained
in the
>>>>         // fcst levels.  If not, print a warning message.
>>>>         if(gc_pd[i].fcst_gci.lvl_type == PresLevel&&
>>>>            gc_pd[i].fcst_gci.lvl_1    !=
gc_pd[i].fcst_gci.lvl_2&&
>>>>            (gc_pd[i].obs_gci.lvl_1<   gc_pd[i].fcst_gci.lvl_1 ||
>>>>             gc_pd[i].obs_gci.lvl_2>   gc_pd[i].fcst_gci.lvl_2)) {
>>>>
>>>>            cout<<  "***WARNING***:
PointStatConfInfo::process_config() ->  "
>>>> <<  "The range of requested observation pressure levels"
>>>> <<  "is not contained within the range of requested"
>>>> <<  "forecast pressure levels.  No vertical interpolation"
>>>> <<  "will be performed for observations falling outside"
>>>> <<  "the range of forecast levels.  Instead, they will be"
>>>> <<  "matched to the single nearest forecast level.\n\n"
>>>> <<  flush;
>>>>
>>>> /I don't see exactly where the calculation is being done in the
code to
>>>> prove or disprove this.  When I ran a comparison between the two
>>>> versions of MET using the same obs and forecast files, I used
ADPSFC
>>>> observations and TMP/P1000-700, which is what you needed to do in
the
>>>> earlier version to make sure you weren't getting rid of surface
>>>> observations at elevation.  With METv0.9, there were 516
observations
>>>> used within my domain of interest, compared to 181 observations
with
>>>> METv3.0.  So it seems like something else is going on that isn't
>>>> immediately clear.  Ideally, I would just use the TMP/Z2, 2-m
diagnosed
>>>> forecast field to compare to the ADPSFC observations (this method
yields
>>>> 690 observations) instead of the P1000-700 method, but not all
the model
>>>> products I use have that diagnosed field.  Anyway....if you could
>>>> clarify what is actually being done with the vertical
interpolation, it
>>>> would definitely put my mind at ease.
>>>>
>>>> Thank you again for all your help.
>>>>
>>>>
>>>> - Jeff Zielonka
>

------------------------------------------------


More information about the Met_help mailing list