Derive and portray Relative Humidity from RTMA data

[alexander-petkov]

This is a subtask for #10.

It took me so long, that now deserves its own ticket.

The task presented some unexpected results (and challenges as well!!), which I will be documenting here.

[alexander-petkov]

From Matt's comment, Relative Humidity is calculated as follows:

RH = (VP(Dewpoint_temperature) / VP(Temperature)) * 100.0

Where:

float CalcVP(int iTemp){ / iTemp: Temperature in Fahrenheit / / Purpose: Calculate the saturation vapor pressure / / Convert Temperature to Kelvin / float kTemp = (iTemp - 32) / 1.8 + 273.15; / Return the saturation vapor pressure at the given temperature / return exp((float)1.81 + (kTemp * 17.27 - 4717.31) / (kTemp - 35.86)); }

[alexander-petkov]

Since the Dewpoint temperature, and Temperature are in Degrees Kelvin in the RTMA files, I didn't need to do conversion to Kelvin, and used cdo with the formula above in an expression, to calculate RH:

./cdo -expr,'2r=(exp(1.81+(2d*17.27- 4717.31) / (2d - 35.86))/exp(1.81+(2t*17.27- 4717.31) / (2t - 35.86)))*100' rtma2p5.t07z.2dvaranl_ndfd.grb2_wexp rhm.grb2

where 2d is Dewpoint Temperature (K), 2t is Temperature (K), and 2r is the variable code for Relative humidity.

[alexander-petkov]

The issue with using cdo

The LCC projection is changed in the resulting RH grib files... I don't know why, especially given the fact that this is a point by point computation.

Original CRS:

 PROJCS["lambert_conformal_conic_1SP", 
  GEOGCS["unknown", 
    DATUM["unknown", 
      SPHEROID["unknown", 6371200.0, 0.0]], 
    PRIMEM["Greenwich", 0.0], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH]], 
  PROJECTION["Lambert_Conformal_Conic_1SP"], 
  PARAMETER["central_meridian", -95.0], 
  PARAMETER["latitude_of_origin", 25.0], 
  PARAMETER["scale_factor", 1.0], 
  PARAMETER["false_easting", 0.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH]]

RH Output CRS:

PROJCS["lambert_conformal_conic_1SP", 
  GEOGCS["unknown", 
    DATUM["unknown", 
      SPHEROID["unknown", 6367470.0, 0.0]], 
    PRIMEM["Greenwich", 0.0], 
    UNIT["degree", 0.017453292519943295], 
    AXIS["Geodetic longitude", EAST], 
    AXIS["Geodetic latitude", NORTH]], 
  PROJECTION["Lambert_Conformal_Conic_1SP"], 
  PARAMETER["central_meridian", -95.0], 
  PARAMETER["latitude_of_origin", 25.0], 
  PARAMETER["scale_factor", 1.0], 
  PARAMETER["false_easting", 0.0], 
  PARAMETER["false_northing", 0.0], 
  UNIT["m", 1.0], 
  AXIS["Easting", EAST], 
  AXIS["Northing", NORTH]]

Essentially, a shifted datum.

[alexander-petkov]

This discrepancy didn't allow me to add RH grib files as an ImageMosaic to Geoserver, altough individual files would display just fine. This had me puzzled for a long time, and esentially took some time debugging in Eclipse to see what is going on.

I defined the new projection in Geoserver's epsg.properties, but I was using gdalinfo to extract CRS information.

GDAL reported it a bit differently, as a LCC with "2 standart parallels":

PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6367470,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",25],
    PARAMETER["standard_parallel_2",25],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",265],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]

This caused a mismatch during ImageMosaic's granule examination, that it caused it to fail, or to index a single granule.

[alexander-petkov]

Debugging in Eclipse showed the actual CRS information, as it was read during configuration, and I defined it with an arbitrary EPSG code:

PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Sphere",6367470,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",25],
    PARAMETER["standard_parallel_2",25],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",265],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]

Also, I had to specify the CRS code for all granules in indexer.xml: <parameter name="MosaicCRS" value="EPSG:45559"/>

Only then all granules were indexed.

[alexander-petkov]

Using GDAL for deriving RH from RTMA Grib files.

It is possible to use GDAL to derive Rel Humidity. I started looking into using GDAL for 2 reasons:

1. The CRS mismatch when using cdo, and the troubles I had configuring these Grib files in Geoserver.
2. It was mentioned that Geotiff might be better suited than Grib, with use of tiling, etc.

Strangely, GDAL reports that Dewpoint Temperature and Air Temperature are in Celsius?

~$ gdalinfo rtma2p5.t02z.2dvaranl_ndfd.grb2_wexp

Metadata:
    GRIB_COMMENT=Temperature [C]

.....

Metadata:
    GRIB_COMMENT=Temperature [C]

So, conversion to Deg in Kelvin is needed when using the RH formula:

gdal_calc.py -of Gtiff -A rtma2p5.t02z.2dvaranl_ndfd.grb2_wexp  --A_band=1 \
    -B rtma2p5.t02z.2dvaranl_ndfd.grb2_wexp --B_band=1 \
    --calc='(exp(1.81+((A+273.15)*17.27- 4717.31) / ((A+273.15) - 35.86))/exp(1.81+((B+273.15)*17.27- 4717.31) / ((B+273.15) - 35.86)))*100' \
    --outfile=rhm.tiff

Edit: avoid converting to Kelvin by setting GRIB_NORMALIZE_UNITS=no

[alexander-petkov]

Current status

The derived RTMA RH data is configured as a mosaic, and sync'ed to the online RTMA archive. I am using cdo to create Grib files for now. This process is now a part of the script which updates the RTMA archive

I have also configured a Mosaic with Geotiff files, using GDAL to make them. These files are uncompressed (faster), untiled, and without adding internal overviews.

I have a dataset with added internal overviews, but haven't succeeded in configuring an ImageMosaic with them, so far... Once I do that, I should run timing metrics for Multidimensional coverages using Grib, Geotiff, and Geotiff with internal overviews.

I think using external (horizontal) tiling, and then adding overviews will be most beneficial. Depending on timing results, I might need to restructure and optimize all weather data for display/querying speed.

To do:

• [ ] Generate uncompressed Geotiffs with RH data. Use a different data type for values b/n 0 and 100. Configure a Mosaic, similar to the one with Grib files.
• [ ] Tile the above Geotiffs, and add internal overviews. Configure an ImageMosaic using this data
• [ ] Collect speed metrics, to measure performance for all 3 types of Mosaics.

Helpful links regarding Raster optimization, and Stress testing with Apache JMeter:

https://geoserver.geo-solutions.it/edu/en/enterprise/raster.html

https://geoserver.geo-solutions.it/edu/en/enterprise/jmeter.html

[alexander-petkov]

Switch to using gdal_calc for deriving RH, as cdo gives unexpected results.