gyanz / pydsstools

Python library for simple HEC-DSS functions
MIT License
80 stars 36 forks source link

Example - writing a GIS Tiff file to DSS. #54

Open LukeWebbJacobs opened 10 months ago

LukeWebbJacobs commented 10 months ago

Hello, I have struggled to write a valid GIS tiff, to a dss file.

I had to leave disabled the 'raise_profile_error', because the values calculated by the pydsstools, appear to be calculated using the top left of the grid, instead of the bottom left.

The comment on the 25th of november here mentions this as well: https://github.com/gyanz/pydsstools/issues/31

The end result is that the raster, ends up shifted up, by the distance of the height of the raster If i provide values that pass your validation.

LukeWebbJacobs commented 10 months ago

Is there a specific example how this should be done?

LukeWebbJacobs commented 10 months ago

Working Example - Conversion from TIFF to DSS Here is an mostly working example to convert a tiff to dss. Similar to your current example 10, writing spatial grid record.

But I have since noticed small projection inconsistancies so any specific example / help would be appreciated!.

Notice this line of code, Which i believe mostly works around the bug. (and causes the profile error to be raised) : ('opt_lower_left_y', ll_y - data.shape[0]),

    from pydsstools.heclib.dss.HecDss import Open
    from pydsstools.heclib.utils import gridInfo
    from pydsstools.core import (lower_left_xy_from_transform)
    import rasterio

    in_tiff = 'in_file.tif'
    out_dss_file = 'a_file.dss'

    with (Open(out_dss_file) as fid):

        # Build output path within dss
        pathname_out = 'a/b/c/d/e/f'

        # Load the raster
        with rasterio.open(in_tiff, 'r') as in_raster:
            # Get data as numpy array
            data = in_raster.read(1)

            # Create dss metadata dictionary
            ll_x, ll_y = lower_left_xy_from_transform(in_raster.transform, data.shape)
            crs_name = in_raster.crs.to_wkt().split('"')[1]

            grid_info = gridInfo()
            grid_info.update([('grid_type', 'specified-time'),
                              ('grid_crs', in_raster.crs.to_wkt()),
                              ('opt_crs_name', crs_name),
                              ('grid_transform', in_raster.transform),
                              ('data_type', 'per-cum'),
                              ('data_units', 'MM'),
                              ('opt_time_stamped', True),
                              ('opt_is_interval', True),

                              ('opt_lower_left_x', ll_x),
                              ('opt_lower_left_y', ll_y - data.shape[0]),
                              ],
                             )

            # Add file to output dss
            fid.put_grid(pathname_out, data, grid_info)