Example - writing a GIS Tiff file to DSS.

[LukeWebbJacobs]

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]

Is there a specific example how this should be done?

[LukeWebbJacobs]

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)

[KmBase]

in grid.pyx , the function lower_left_xy_from_transform is :

def lower_left_xy_from_transform(transform,shape,cell_zero_xcoord=0,cell_zero_ycoord=0):
    cdef:
        float xmin,ymax
        float ymin # verbose to make calculation more clear
        int lower_left_x,lower_left_y
    cellsize = transform[0]
    cellsize_y = transform[4]
    rows,cols = shape
    xmin = transform[2] - cell_zero_xcoord
    ymax = transform[5] - cell_zero_ycoord
    ymin = ymax + rows * cellsize_y
    lower_left_x = <int>(floor(xmin/cellsize))
    lower_left_y = <int>(floor(ymax/cellsize))
    return (lower_left_x,lower_left_y)

why the value of lower_left_y is y_max? There is a source about the coordinate-reference-systems

For identification, each cell in the grid has a pair of integer indices (i, j) indicating the position, by cell count, of its southwest (or minimum-x, minimum-y) corner, relative to the UTM zone coordinate origin. To find the indices of the cell in which a point is located, find the point’s easting and northing in the projected coordinate system defined above, and calculates the indices with the following formulas. i = floor( easting / cellsize ) j = floor( northing / cellsize ) where floor(x) is the largest integer less than or equal to x.

We know that transform[2] and transform[5] is the upper_left location. So the xmin and ymin are absolutely correct in the lower_left_xy_from_transform function.But why ymin is unused.

[KmBase]

By the way, I find that when x_cellsize is not equal to y_cellsize, the result will be far away from the result generated by votex.

[ReasonableRadish]

I'm finding this issue still with the most recent version of Pydsstools. It seems even the example 10 given on the front page has this same issue with the lower_left_xy_from_transform such that the computed lower left coords aren't matching the dataset.

Output from example 10: "ERROR:root:opt_lower_left_x, opt_lower_left_y has issue or both are incorrect Given = (0, 0), computed = (10, 50)".

Related to KmBase's comment about cellsize, I've found that with the affine transform the x and y cell sizes are the same but opposite signs (transform[4] is negative) I assume this is a matrix algebra thing to put it simply

[LukeWebbJacobs]

It would be great for someone to solve this. I use pydss everyday, and since raising the issue, have now created over 1 million dss files, all with the grids in slightly the wrong location :)

[ReasonableRadish]

Hm your grids still write with that error? I found that I was able to write the grids to the dss and it has the meta data and what not but doesn't seem to transfer the values correctly because its showing up as all NaNs. I can click to view the dss entry but nothing appears in the grid. Can I ask what you're writing all these dss files for? Are you writing from a tif?

[LukeWebbJacobs]

I download .GRB (Gridded binary data) weather models. And convert them to dss to do some HEC modelling multiple times a day.

However, the above tif code, should work for you.

The issue, is if I use dss, and convert back to tiff I can see all my grids have slightly shifted.

I am unable to create a unit test that:

• Load tiff file
• Save to Dss
• Load Dss
• Save to Tiff
• Check input tiff == output tiff.

I am of course, able to do this test, for my GRB files, and tiff files :)

Instead, I know the scale at which my tiffs are created, so know roughly how misalinged the data will be to my truth, and just add that to my error tolerances :D (And also buffer slightly to ensure all my extent has grid data after shitfing)

[ReasonableRadish]

Oh ok, interesting thanks for sharing! If you don't mind, I would be interested to see an example if you had one to share which shows the grib files you've used and snippets of the functioning code. I'm still trying to get my tif code to work using the CA/NV regional QPF data. I seem to almost have it working, but it's still writing an empty dss file.

[LukeWebbJacobs]

Its quite bespoke code, here is an example. I cant include any grib data sorry as it is proprietry, but at this point in the code, the data is the same as any old raster.

[LukeWebbJacobs]

def add_forecast_to_dss(clipped_gri.txt

[gyanz]

@LukeWebbJacobs Would revising the computation of lower_left_y in lower_left_xy_from_transform address this issue? For cell size, DSS stores a single value for Specified grid, implying square cells, so I assumed equal x and y cell sizes.

[LukeWebbJacobs]

@gyanz Yes I think it will help.

Maybe my cell sizes were not equal or maybe something else I cant find in an hour or 2 debugging!, I stopped working on this now.

But using version 4.2 I think my transforms seem ok now and my raster converts both ways. (Did you fix this already?)

[gyanz]

@LukeWebbJacobs Yes I have fixed it. I have made further fixes in the 3.0.0 beta version.