HydrologicEngineeringCenter / Vortex

data processing utilities
MIT License
25 stars 7 forks source link

Part D and Part E fields not generating properly in output DSS file #91

Closed jlgutenson closed 1 year ago

jlgutenson commented 1 year ago

Hello,

I've been attempting to use the Vortex API and Jython to convert gauge corrected MRMS data into a clipped DSS file to import in HEC-HMS.

After trying to use the resulting DSS file in my simulation with no luck, I compared the resulting DSS file to one generated through the HEC-HMS interface that was simulating properly.

Upon comparison in HEC-DSSVue, I noticed that Part D and Part E fields in the DSS file are not being generated properly. Here's my script below, any suggestions?

from mil.army.usace.hec.vortex.io import BatchImporter
from mil.army.usace.hec.vortex.geo import WktFactory
import os

met_grib_dir = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/Scripts/build_hms_inputs/mrms_gage_corrected"

clip_shp = "C:/Users/Joseph Gutenson/Desktop/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1. RVD/Extrass/RVDWshed.shp"

destination = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1_RVD/1_HEC-HMS_Model/RVD_NAD8310171/RVDJune2018_JLG_scripted.dss"

# see if dss file already exists and if so, remove it
if os.path.exists(destination):
    os.remove(destination)
else:
    pass

# list all grib files in the directory 
in_files = os.listdir(met_grib_dir)

# pair each grib file with its full path and pack into a list
in_files_path = []
for in_file in in_files:
    if in_file.endswith("grib2.gz"):
        in_file_full_path = os.path.join(met_grib_dir,in_file)
        in_files_path.append(in_file_full_path)
    else:
        pass

# sort the list of grib files, doesn't seem to make a difference
in_files_path.sort()

# parameters and set up for Vortex to create the dss file
# documented here: https://github.com/HydrologicEngineeringCenter/Vortex/wiki/Batch-Import-Options
variables = ['GaugeCorrQPE01H_altitude_above_msl']

geo_options = {
    'pathToShp': clip_shp,
    'targetCellSize': '2000',
    'targetEpsg': '26914',
    'resamplingMethod': 'Bilinear'
}

write_options = {
    'partA': 'SHG',
    'partB': 'RVD',
    'partC': 'PRECIPITATION',
    'partF': 'MRMS',
    'dataType': 'PER-CUM',
    'units': 'MM'
}

myImport = BatchImporter.builder() \
                .inFiles(in_files_path) \
                .variables(variables) \
                .geoOptions(geo_options) \
                .destination(destination) \
                .writeOptions(write_options) \
                .build()

myImport.process()

Thanks, Joseph

jlgutenson commented 1 year ago

image Here's the resulting DSS in DSSVue.

tombrauer commented 1 year ago

Hi @jlgutenson, I see that the start end time are the same. Not all datasets are good about declaring the duration over which the data applies and if I recall correctly MRMS is one of those datasets. There is logic in Vortex that tries to parse the name of the file and if it is MRMS, apply the appropriate temporal extent. This is not ideal, but it works. There must be something about your case where this breaks down. As a quick fix you can use the time-shifter to shift just the start dates or end dates by one hour. The data will then become period data and that should make HEC-HMS happy. If you can zip and post a few of the files you are writing your script against, I can take a look and try to figure out what is happening.

jlgutenson commented 1 year ago

Thanks @tombrauer! Here's a zip of a few of those files. I'll see if I can play with things on my end. mrms_gage_corrected.zip

jlgutenson commented 1 year ago

What's strange is that if I go through the HEC-HMS GUI import grid process with the same files, the output DSS generates correctly.

jlgutenson commented 1 year ago

I seemed to have resolved the issue with Part D and Part E being out of order. I decompressed the .gz file and then used the resulting .grib2 file.

However, the output DSS file is still roughly half of the size of the DSS file created through the HEC-HMS GUI and HEC-HMS still doesn't like it (HEC-HMS just aborts the run).

I also notice that I'm missing one time step in my Jython/Vortex-API DSS file that is in the DSS file created by the HEC-HMS GUI.

jlgutenson commented 1 year ago

Current code for reference, using Vortex 0.11.0:

# third party Java modules
from mil.army.usace.hec.vortex.io import BatchImporter
from mil.army.usace.hec.vortex.geo import WktFactory
from mil.army.usace.hec.vortex.io import DataReader
from mil.army.usace.hec.vortex.math import Shifter
from java.time import Duration

# native modules to add
import os
import shutil
import gzip

# list of variables that we'll eventually convert to inputs in a function.
met_grib_dir = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/Scripts/build_hms_inputs/mrms_gage_corrected"
clip_shp = "C:/Users/Joseph Gutenson/Desktop/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1. RVD/Extrass/RVDWshed.shp"
destination_1 = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1_RVD/1_HEC-HMS_Model/RVD_NAD8310171/RVDJune2018_JLG_scripted_1.dss"
# destination_2 = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1_RVD/1_HEC-HMS_Model/RVD_NAD8310171/RVDJune2018_JLG_scripted_2.dss"

# see if dss file already exists and if so, remove it
if os.path.exists(destination_1):
    os.remove(destination_1)
else:
    pass

# list all grib files in the directory 
in_files = os.listdir(met_grib_dir)

# pair each grib file with its full path and pack into a list
in_files_path = []
for in_file in in_files:
    if in_file.endswith("grib2.gz"):
        # decompress *.gz files
        in_file_full_path = os.path.join(met_grib_dir,in_file)
        with gzip.open(in_file_full_path, 'rb') as f_in:
            out_file_full_path = in_file_full_path[:-3]
            with open(out_file_full_path, 'wb') as f_out:
                shutil.copyfileobj(f_in, f_out)
        in_files_path.append(out_file_full_path)
    else:
        pass

# sort the list of grib files, doesn't seem to make a difference
in_files_path.sort()

# parameters and set up for Vortex to create the dss file
# documented here: https://github.com/HydrologicEngineeringCenter/Vortex/wiki/Batch-Import-Options
variables = ['GaugeCorrQPE01H_altitude_above_msl']

geo_options = {
    'pathToShp': clip_shp,
    'targetCellSize': '2000',
    'targetEpsg': '26914',
    'resamplingMethod': 'Bilinear'
}

write_options = {
    'partA': 'SHG',
    'partB': 'RVD',
    'partF': 'MRMS',
    'dataType': 'PER-CUM',
    'units': 'MM'
}

myImport = BatchImporter.builder() \
                .inFiles(in_files_path) \
                .variables(variables) \
                .geoOptions(geo_options) \
                .destination(destination_1) \
                .writeOptions(write_options) \
                .build()

myImport.process()
tombrauer commented 1 year ago

I'm able to run your script with the sample files and generate seemingly valid records in an HEC-DSS file. Are you able to determine what is wrong with the script-generated file other than the missing record? It may help to delete the existing HEC-DSS file and let your script recreate it. If records are overwritten many times over, that can lead to file size discrepancies.

jlgutenson commented 1 year ago

One other difference I noticed is that DSSVue labels the file version differently. The DSS file that works in HMS is labeled as File Version: 7-IQ. The file version of the Jython generated DSS file is 7-IP. Could that be causing some issues with the HMS simulations?

I've also noticed that the file size of the Jython generated DSS file is about half the size of the HMS generated DSS file.

tombrauer commented 1 year ago

I wouldn't expect the DSS file version to matter. If you want to email a copy of your project with the failing simulation, I can take a look to see why it is failing.

jlgutenson commented 1 year ago

FIgured it out @tombrauer and, as I suspected, it was something simple and silly.

I was specifying the EPSG to match the HMS model's coordinate system. Instead, I passed the 'targetWkt': WktFactory.shg() in the geo_options and that did the trick. I'm not quite sure what that does but it gives the DSS file an Albers coordinate system.

Here's the completed script below for posterity:


from mil.army.usace.hec.vortex.io import BatchImporter
from mil.army.usace.hec.vortex.geo import WktFactory
from mil.army.usace.hec.vortex.io import DataReader
from mil.army.usace.hec.vortex.math import Shifter

# native modules to add
import os
import shutil
import gzip

# list of variables that we'll eventually convert to inputs in a function.
met_grib_dir = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/Scripts/build_hms_inputs/mrms_gage_corrected"
clip_shp = "C:/Users/Joseph Gutenson/Desktop/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1. RVD/Extrass/RVDWshed.shp"
destination_1 = "D:/Gutenson_RATES/TWDB-FIF-LRGVDC/2023/1.2.2.2.2/Models/HEC_HMS_411beta16/1_RVD/1_HEC-HMS_Model/RVD_NAD8310171/RVDJune2018_JLG_scripted_1.dss"

# see if dss file already exists and if so, remove it
if os.path.exists(destination_1):
    os.remove(destination_1)
else:
    pass

# list all grib files in the directory 
in_files = os.listdir(met_grib_dir)

# pair each grib file with its full path and pack into a list
in_files_path = []
for in_file in in_files:
    if in_file.endswith("grib2.gz"):
        # decompress *.gz files
        in_file_full_path = os.path.join(met_grib_dir,in_file)
        with gzip.open(in_file_full_path, 'rb') as f_in:
            out_file_full_path = in_file_full_path[:-3]
            with open(out_file_full_path, 'wb') as f_out:
                shutil.copyfileobj(f_in, f_out)
        in_files_path.append(out_file_full_path)
    else:
        pass

# sort the list of grib files, doesn't seem to make a difference
in_files_path.sort()

# parameters and set up for Vortex to create the dss file
# documented here: https://github.com/HydrologicEngineeringCenter/Vortex/wiki/Batch-Import-Options
variables = ['GaugeCorrQPE01H_altitude_above_msl']

geo_options = {
    'pathToShp': clip_shp,
    'targetCellSize': '2000',
    'targetWkt': WktFactory.shg(),
    'resamplingMethod': 'Bilinear'
}

write_options = {
    'partA': 'SHG',
    'partB': 'RVD',
    'partF': 'MRMS',
    'dataType': 'PER-CUM',
    'units': 'MM'
}

myImport = BatchImporter.builder() \
                .inFiles(in_files_path) \
                .variables(variables) \
                .geoOptions(geo_options) \
                .destination(destination_1) \
                .writeOptions(write_options) \
                .build()

myImport.process()
tombrauer commented 1 year ago

Good to know! WktFactory.shg() will simply return a WKT representation for the Albers Equal Area Conic/SHG projection, a.k.a. EPSG: 5070.

    public static String getShg() {
        return "PROJCS[\"USA_Contiguous_Albers_Equal_Area_Conic_USGS_version\"," +
                "        GEOGCS[\"GCS_North_American_1983\"," +
                "        DATUM[\"D_North_American_1983\"," +
                "        SPHEROID[\"GRS_1980\",6378137.0,298.257222101]]," +
                "PRIMEM[\"Greenwich\",0.0]," +
                "UNIT[\"Degree\",0.0174532925199433]]," +
                "PROJECTION[\"Albers\"]," +
                "        PARAMETER[\"False_Easting\",0.0]," +
                "PARAMETER[\"False_Northing\",0.0]," +
                "PARAMETER[\"Central_Meridian\",-96.0]," +
                "PARAMETER[\"Standard_Parallel_1\",29.5]," +
                "PARAMETER[\"Standard_Parallel_2\",45.5]," +
                "PARAMETER[\"Latitude_Of_Origin\",23.0]," +
                "UNIT[\"Meter\",1.0]]";
    }

So your code was previously reprojecting to UTM 14 N and is now reprojecting to SHG.

jlgutenson commented 1 year ago

Hey @tombrauer, does it make a difference if the grib2 data you're importing is in a projected coordinated system (e.g., like HRRR Grib2 data)?

tombrauer commented 1 year ago

@jlgutenson that should be fine. I would just be sure to reproject to a projected coordinate system as part of the import process.

jlgutenson commented 1 year ago

So if I'm importing a grib with a Lambert Conformal Conic projection, with the horizontal units in kilometers, the BatchImporter handles that behind the scenes?

tombrauer commented 1 year ago

It should... there's only one way to find out 😄

jlgutenson commented 1 year ago

I gave the ole' sporting try. The import ran suspiciously fast and the resulting DSS grid did have any values in DSSVue. That makes me think that Vortex thinks my watershed is somewhere in the Indian Ocean (we all know that result lol). Everything else in the DSS file seems correct though. I'll see if I can't convert the grib from Lambert projection to something geographic (i.e., NAD83).

tombrauer commented 1 year ago

Gotcha, horizontal units in kilometers is not very common, but I thought Vortex could handle it. Feel free to post your data source and I can take a look. Let's do it on another issue. I think this one is exhausted.

jlgutenson commented 1 year ago

Sounds good, I'll close this to get it out of your hair.