search

JCSDA / CRTMv3

CRTMv3 repository for coordinated development and releases. Code history is not carried in this repository prior to v3, to reduce the cloning overhead. For v2.x history leading up to v3, see JCSDA/crtm repository.
Other
6 stars 6 forks source link

Create CRTM_RTSolution_ReadFile module for netCDF format #72

Closed chengdang closed 9 months ago

chengdang commented 1 year ago

Create CRTM_RTSolution_ReadFile module for netCDF format

chengdang commented 1 year ago

Hello @BenjaminTJohnson. Here are some thoughts on this issue. Please let me know your suggestions.

After reviewing the data structure and some initial module implementation, I realize that there will be some relatively substantial changes required if we want to replace all binary i/o in ctests. Mostly because the data structure in binary/netCDF files are slightly different.

In binary output files, the RTSolution are "stacked" based on the number of channels and profiles. In each of these output blocks indexed by channel and profile, the number of layers in one profile, and the number of stokes do not matter. For each ctests, there is only one output file needed for all channels/profiles being tested. The binary output structure aligns well with the RTSolution data structure used in CRTM: RTSolution(n_Channels, n_Profiles) https://github.com/JCSDA/CRTMv3/blob/3ec6c9a96bf55281a941c8af11fc248ee7d8cac0/src/RTSolution/CRTM_RTSolution_Define.f90#L1290

In netCDF output files, I rearranged RTSolution so that for each profile, there is one output NC file containing all channels. This is chosen mostly because the number of layers can be different across the input profiles. For this I/O routines, the number of layers in a profile and the number of stokes need to be included in the file. Such netCDF data structure is more convenient for offline testing and data visualization. https://github.com/JCSDA/CRTMv3/blob/3ec6c9a96bf55281a941c8af11fc248ee7d8cac0/src/RTSolution/CRTM_RTSolution_Define.f90#L1460

It is still possible to replace all binary i/o in ctests, while before which I'd like to double check and make sure the netCDF output structure is what we want. One modification that we can do is to assume the number of layers are the same in all profiles (per test) that allows us to generate one file for each sensor (multiple profiles, all channels) - similar to the Binary output file.

Here is the NC output data for reference:

netcdf atms_npp_Profile1 {
dimensions:
    n_Layers = 92 ;
    n_Channels = 22 ;
    n_Stokes = 1 ;
variables:
    int Sensor_Channel(n_Channels) ;
        Sensor_Channel:units = "N/A" ;
        Sensor_Channel:_FillValue = 0 ;
    int n_Full_Streams(n_Channels) ;
        n_Full_Streams:units = "N/A" ;
        n_Full_Streams:_FillValue = 0 ;
    double SSA_Max(n_Channels) ;
        SSA_Max:units = "fraction (0->1)" ;
        SSA_Max:_FillValue = -999. ;
    double SOD(n_Channels) ;
        SOD:units = "1" ;
        SOD:_FillValue = -999. ;
    double Surface_Emissivity(n_Channels) ;
        Surface_Emissivity:units = "fraction (0->1)" ;
        Surface_Emissivity:_FillValue = -999. ;
    double Surface_Reflectivity(n_Channels) ;
        Surface_Reflectivity:units = "fraction (0->1)" ;
        Surface_Reflectivity:_FillValue = -999. ;
    double Up_Radiance(n_Channels) ;
        Up_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Up_Radiance:_FillValue = -999. ;
    double Down_Radiance(n_Channels) ;
        Down_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Down_Radiance:_FillValue = -999. ;
    double Down_Solar_Radiance(n_Channels) ;
        Down_Solar_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Down_Solar_Radiance:_FillValue = -999. ;
    double Surface_Planck_Radiance(n_Channels) ;
        Surface_Planck_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Surface_Planck_Radiance:_FillValue = -999. ;
    double Total_Cloud_Cover(n_Channels) ;
        Total_Cloud_Cover:units = "fraction (0->1)" ;
        Total_Cloud_Cover:_FillValue = -999. ;
    double R_clear(n_Channels) ;
        R_clear:units = "fraction (0->1)" ;
        R_clear:_FillValue = -999. ;
    double Tb_clear(n_Channels) ;
        Tb_clear:units = "Kelvin" ;
        Tb_clear:_FillValue = -999. ;
    double Radiance(n_Channels) ;
        Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Radiance:_FillValue = -999. ;
    double Brightness_Temperature(n_Channels) ;
        Brightness_Temperature:units = "Kelvin" ;
        Brightness_Temperature:_FillValue = -999. ;
    double Solar_Irradiance(n_Channels) ;
        Solar_Irradiance:units = "Watts per Square Metre per Micron (W m^-2 micron^-1)" ;
        Solar_Irradiance:_FillValue = -999. ;
    double Stokes(n_Stokes, n_Channels) ;
        Stokes:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Stokes:_FillValue = -999. ;
    double Upwelling_Overcast_Radiance(n_Layers, n_Channels) ;
        Upwelling_Overcast_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Upwelling_Overcast_Radiance:_FillValue = -999. ;
    double Upwelling_Radiance(n_Layers, n_Channels) ;
        Upwelling_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Upwelling_Radiance:_FillValue = -999. ;
    double Layer_Optical_Depth(n_Layers, n_Channels) ;
        Layer_Optical_Depth:units = "1" ;
        Layer_Optical_Depth:_FillValue = -999. ;
    double Single_Scatter_Albedo(n_Layers, n_Channels) ;
        Single_Scatter_Albedo:units = "fraction (0->1)" ;
        Single_Scatter_Albedo:_FillValue = -999. ;

// global attributes:
        :Sensor_ID = "atms_npp" ;
}
chengdang commented 1 year ago

Talked with Fabio that we could safely assume the number of layers are the same for all profiles when CRTM is called. Will modify these modules to generate one RTS output file per sensor.

chengdang commented 1 year ago

Initial ctests passed with NetCDF I/O. Need to add more global attributes require for RTS comparison.

chengdang commented 1 year ago

Finished updating/developing all modules. Modified only unit test test_forward to show how to use these modules in CRTM. The RTS files in netCDF format are much larger than those in binary format (up to 10 times bigger). More discussion are needed on if these interfaces should replace the existing binary ones.

As an example, here are the data structures in file test/results/forward/test_Simple_v.abi_gr.RTSolution.nc:

netcdf test_Simple_v.abi_gr.RTSolution {
dimensions:
    n_Profiles = 2 ;
    n_Layers = 92 ;
    n_Channels = 6 ;
    n_Stokes = 1 ;
variables:
    int Sensor_Channel(n_Channels) ;
        Sensor_Channel:units = "N/A" ;
        Sensor_Channel:_FillValue = 0 ;
    int n_Full_Streams(n_Profiles, n_Channels) ;
        n_Full_Streams:units = "N/A" ;
        n_Full_Streams:_FillValue = 0 ;
    double SSA_Max(n_Profiles, n_Channels) ;
        SSA_Max:units = "fraction (0->1)" ;
        SSA_Max:_FillValue = -999. ;
    double SOD(n_Profiles, n_Channels) ;
        SOD:units = "1" ;
        SOD:_FillValue = -999. ;
    double Surface_Emissivity(n_Profiles, n_Channels) ;
        Surface_Emissivity:units = "fraction (0->1)" ;
        Surface_Emissivity:_FillValue = -999. ;
    double Surface_Reflectivity(n_Profiles, n_Channels) ;
        Surface_Reflectivity:units = "fraction (0->1)" ;
        Surface_Reflectivity:_FillValue = -999. ;
    double Up_Radiance(n_Profiles, n_Channels) ;
        Up_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Up_Radiance:_FillValue = -999. ;
    double Down_Radiance(n_Profiles, n_Channels) ;
        Down_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Down_Radiance:_FillValue = -999. ;
    double Down_Solar_Radiance(n_Profiles, n_Channels) ;
        Down_Solar_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Down_Solar_Radiance:_FillValue = -999. ;
    double Surface_Planck_Radiance(n_Profiles, n_Channels) ;
        Surface_Planck_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Surface_Planck_Radiance:_FillValue = -999. ;
    double Total_Cloud_Cover(n_Profiles, n_Channels) ;
        Total_Cloud_Cover:units = "fraction (0->1)" ;
        Total_Cloud_Cover:_FillValue = -999. ;
    double R_clear(n_Profiles, n_Channels) ;
        R_clear:units = "fraction (0->1)" ;
        R_clear:_FillValue = -999. ;
    double Tb_clear(n_Profiles, n_Channels) ;
        Tb_clear:units = "Kelvin" ;
        Tb_clear:_FillValue = -999. ;
    double Radiance(n_Profiles, n_Channels) ;
        Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Radiance:_FillValue = -999. ;
    double Brightness_Temperature(n_Profiles, n_Channels) ;
        Brightness_Temperature:units = "Kelvin" ;
        Brightness_Temperature:_FillValue = -999. ;
    double Solar_Irradiance(n_Profiles, n_Channels) ;
        Solar_Irradiance:units = "Watts per Square Metre per Micron (W m^-2 micron^-1)" ;
        Solar_Irradiance:_FillValue = -999. ;
    double Stokes(n_Profiles, n_Stokes, n_Channels) ;
        Stokes:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Stokes:_FillValue = -999. ;
    double Upwelling_Overcast_Radiance(n_Profiles, n_Layers, n_Channels) ;
        Upwelling_Overcast_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Upwelling_Overcast_Radiance:_FillValue = -999. ;
    double Upwelling_Radiance(n_Profiles, n_Layers, n_Channels) ;
        Upwelling_Radiance:units = "Watts per Square Metre per Micron per Rad (W m^-2 micron^-1 rad^-1)" ;
        Upwelling_Radiance:_FillValue = -999. ;
    double Layer_Optical_Depth(n_Profiles, n_Layers, n_Channels) ;
        Layer_Optical_Depth:units = "1" ;
        Layer_Optical_Depth:_FillValue = -999. ;
    double Single_Scatter_Albedo(n_Profiles, n_Layers, n_Channels) ;
        Single_Scatter_Albedo:units = "fraction (0->1)" ;
        Single_Scatter_Albedo:_FillValue = -999. ;

// global attributes:
        :Sensor_ID = "v.abi_gr" ;
        :WMO_Satellite_ID = 1023 ;
        :WMO_Sensor_ID = 2047 ;
        :RT_Algorithm_Name = "ADA" ;
}
chengdang commented 1 year ago

Code updates have been pushed to branch https://github.com/JCSDA/CRTMv3/tree/feature/cd_RTS_netCDF_IO