search

alexander-petkov / wfas

A placeholder for the WFAS project.
4 stars 1 forks source link

Add Climate Forecast System v2 Data #50

Open alexander-petkov opened 2 years ago

alexander-petkov commented 2 years ago

Explore adding the Climate Forecast System V2 Forecast into GeoServer Same variables as before https://cfs.ncep.noaa.gov/cfsv2/downloads.html https://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/

alexander-petkov commented 2 years ago

A detailed PowerPoint presentation with the contents of each dataset https://cfs.ncep.noaa.gov/cfsv2.info/Operational.CFSv2.info.ppt

alexander-petkov commented 2 years ago

Variables, and their datasets

Full Name Abreviation Dataset Units Band # (in grb2 file)
Temperature 2t flx deg K 38
Relative Humidity 2r pgb % 368
Total Cloud cover tcc flx % 53
Downward Short-Wave Rad. Flux dswrf flx W/(m^2) 16
Wind:
1. u-component of wind
2. v-component of wind
10u
10v		 flx [m/s]
36
37
Precipitation rate prate flx [kg/m^2/s] 31
alexander-petkov commented 2 years ago

We need data up to 2 weeks out, as with GFS Files are 6-hourly, containing a single timestep.

alexander-petkov commented 2 years ago

Getting the PGB list of files for a 2-week span:

curl -s --list-only https://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod/cfs/cfs.20211116/00/6hrly_grib_01/ |grep  -oP '(?<=href=")pgbf.*.grb2(?=")'|sort|head -n 56

The forecast is for Nov 16,2021

alexander-petkov commented 2 years ago

Consider calculating Relative Humidity

Consider calculating RH from Specific Humidity and Surface Pressure from the flx dataset. Relative humidity at 2m above ground is found in the pgb dataset, but each file is 21MB. This seems costly with regards to bandwidth, just to extract one variable.

In adition to Temperature at 2m, the FLX dataset also contains Specific humidity and Surface Pressure. The same calculation (as with the GMAO dataset) can be used to derive RH:

gdal_calc.py --format=GTiff -A input.tiff --A_band=1 \
   -B input.tiff --B_band=2 \ 
   -C input.tiff --C_band=3 \
   --calc='( C * B / (0.378 * C + 0.622))/(6.112 * exp((17.67 * (A-273.15))/(A-29.65)))' \
   --outfile=rh.tif

Temperature at 2m: Band 38 Specific Humidity at 2m: Band 39 Pressure at ground or water surface: Band 40

Reference: R metutils

Results can then be compared with 2r in the PGB dataset.

Edit: Revisit this:

Results can then be compared with 2r in the PGB dataset.

alexander-petkov commented 2 years ago

Original source data is referenced as 0-360 deg. Rewrite to -180 to 180 deg

alexander-petkov commented 2 years ago

Re-writing data to -180 to 180 grid with cdo looses the 10v field.

The grb2 files are peculiar--cdo doesn't list the 10v band, but gdal does.

I might have to write a custom routine to swap hemispheres and center data around the Greenwich meridian.

alexander-petkov commented 2 years ago

I crafted a VRT, which defines a Geotransform starting from ~180, and also swaps hemispheres:

<VRTDataset rasterXSize="384" rasterYSize="190">
  <SRS dataAxisToSRSAxisMapping="2,1">GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]</SRS>
  <GeoTransform> -1.8046849864726000e2,  9.3749934713541672e-01,  0.0000000000000000e+00,  8.9750684199999995e+01,  0.0000000000000000e+00, -9.4736675894736833e-01</GeoTransform>
  <Metadata>
    <MDI key="AREA_OR_POINT">Area</MDI>
  </Metadata>
  <VRTRasterBand dataType="Int16" band="1">
    <Metadata>
      <MDI key="STATISTICS_MAXIMUM">100</MDI>
      <MDI key="STATISTICS_MEAN">78.142420504386</MDI>
      <MDI key="STATISTICS_MINIMUM">0</MDI>
      <MDI key="STATISTICS_STDDEV">12.646907077023</MDI>
      <MDI key="STATISTICS_VALID_PERCENT">100</MDI>
    </Metadata>
    <NoDataValue>-9999</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">tif/202111220000.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="384" RasterYSize="190" DataType="Int16" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="192" ySize="190" />
      <DstRect xOff="192" yOff="0" xSize="192" ySize="190" />
    </SimpleSource>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">tif/202111220000.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="384" RasterYSize="190" DataType="Int16" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="192" yOff="0" xSize="192" ySize="190" />
      <DstRect xOff="0" yOff="0" xSize="192" ySize="190" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

The upper x for the Geotransform was calculated by using the upper x of the original data as follows:

>>> -0.4687493- (384/2)*0.937498694516971
-180.46849864725846

Here is an illustration of the initial problem (grayscale), and desired result portrayed by the VRT (pseudocolor): Screenshot from 2021-11-30 23-43-53

The initial 0-360 deg grayscale raster is overlayed over the pseudocolored VRT.

Sources: https://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial https://web.archive.org/web/20130618003511/http:/www1.eonfusion.com/manual/index.php/Manipulate_rasters_with_GDAL_VRT

wmjolly commented 2 years ago

Very nice solution Alex.

On Tue, Nov 30, 2021, 23:57 alexander-petkov @.***> wrote:

I crafted a VRT, which defines a Geotransform starting from ~180, and also swaps hemispheres:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]] -1.8046849864726000e2, 9.3749934713541672e-01, 0.0000000000000000e+00, 8.9750684199999995e+01, 0.0000000000000000e+00, -9.4736675894736833e-01 Area 100 78.142420504386 0 12.646907077023 100 -9999 Gray tif/202111220000.tif 1 tif/202111220000.tif 1

Here is an illustration of the initial problem (grayscale), and desired result portrayed by the VRT (pseudocolor): [image: Screenshot from 2021-11-30 23-43-53] https://user-images.githubusercontent.com/39599557/144185223-d1737eae-a013-42eb-8825-bfeeaa8516d1.png

The initial 0-360 deg grayscale raster is overlayed over the pseudocolored VRT.

Sources: https://trac.osgeo.org/gdal/wiki/UserDocs/RasterProcTutorial

https://web.archive.org/web/20130618003511/http:/www1.eonfusion.com/manual/index.php/Manipulate_rasters_with_GDAL_VRT

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/alexander-petkov/wfas/issues/50#issuecomment-983344394, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4G3D2FOJ2AO7RV3SZ4QMTUOXBMZANCNFSM5IEA3NZA .

alexander-petkov commented 2 years ago

The 7 weather variables from CFS data are now added to AWS Geoserver. Screenshot 2021-12-02 12 15 05

alexander-petkov commented 2 years ago

Comitted update script for CFS data.

alexander-petkov commented 2 years ago

I reverted the CFS update script to previous version for EBS storage instead of S3.

alexander-petkov commented 4 months ago

Added derive_prate() function to transform Total precip from [kg/m^2 s] to [kg/m^2] in commit 79c2fcd

@wmjolly Question: Isn't it better to leave the original units, as 6-hr accumulation is perhaps not very intuitive?