WPSDataset: Wraps GDAL dataset and provides common wps operations

[conbrad]

Allows us to:

• multiply rasters with the existing multiply function: output_ds = a_ds * b_ds
• use with to open rasters within the context and autoclose them outside the context to prevent memory leaks:

  with WPSDataset('path/or/key/to/dataset') as wps_ds:
      # do stuff with wps_ds
  
  # wps_ds underlying gdal dataset is now closed
  

Test Links:

Landing Page
MoreCast
Percentile Calculator
C-Haines
FireBat
FireBat bookmark
Auto Spatial Advisory (ASA)
HFI Calculator

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 80.56%. Comparing base (9466769) to head (c982144).
Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #4010      +/-   ##
==========================================
- Coverage   80.60%   80.56%   -0.04%     
==========================================
  Files         302      303       +1     
  Lines       11456    11564     +108     
  Branches      549      549              
==========================================
+ Hits         9234     9317      +83     
- Misses       2104     2129      +25     
  Partials      118      118              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[brettedw]

Looking good, I like the idea!

[brettedw]

Looks great, good idea!

[dgboss]

Maybe @brettedw can chime in, but I think there are some subtleties we need to be cautious of here when setting the extent and xRes/yRes in WarpOptions.

If we're warping a raster with a geographic coordinate system (think coordinates are lat/lon in degrees) to a projected coordinate system (coordinates are in feet or meters as measured from an arbitrary point) and vice versa, we're going to run into problems I think.

[conbrad]

Interesting, this was pulled over from

wps/api/app/utils/geospatial.py

Line 43 in 6b15aa7

return gdal.Warp(output_path, source_ds, dstSRS=ds_to_match.GetProjection(), outputBounds=extent, xRes=x_res, yRes=y_res, resampleAlg=resample_method.value)

that we currently use.

[dgboss]

Yeah, I didn't notice the issue with the original implementation. I caught it as I was generating those test files. I created the first test file in EPSG4326 which uses lat/lon coordinates and has x and y resolution measured in degrees. I basically used the warp_to_match code to then create the EPSG3857 test file but I kept getting errors about not having enough disk space. The problem was that I was taking the x/y res in degrees which was very small (0.003333333333067) and passing that in to gdal.Warp which then tries to create an output raster where the x/y resolution is 0.003333333333067 metres!
Additionally, I ran into problems with the 4326 extent not making sense in the 3857 spatial reference.

[dgboss]

I'm going to do a quick test and see if specifying a source spatial reference works around the issue.

[dgboss]

Boo...specifying the srcSRS WarpOption doesn't help.
ERROR 3: 3005_lats.tif: Free disk space available is 478468026368 bytes, whereas 635088851384008 are at least necessary.

[brettedw]

Riiiiiight, ya that's an interesting one. Gonna have to think about a solution to that one. Maybe for now we can just implement some sort of check to see if both datasets are in geographic or projected. My brain can't do that right now but maybe we have something like this to use as a check (I don't know if this works, came from chatgippity):

def is_geographic(crs_wkt: str) -> bool:
    """
    Check if the given CRS (WKT string) is geographic (i.e., using degrees).
    :param crs_wkt: The WKT string of the CRS.
    :return: True if geographic, False if projected.
    """
    srs = osr.SpatialReference()
    srs.ImportFromWkt(crs_wkt)
    return srs.IsGeographic()

def get_crs_units(crs_wkt: str) -> str:
    """
    Get the units of the CRS (either 'degrees' for geographic, or 'meters/feet' for projected).
    :param crs_wkt: The WKT string of the CRS.
    :return: The units of the CRS ('degrees', 'meters', 'feet', etc.)
    """
    srs = osr.SpatialReference()
    srs.ImportFromWkt(crs_wkt)
    
    if srs.IsGeographic():
        return 'degrees'
    else:
        return srs.GetLinearUnitsName()

Otherwise maybe we have to do some units conversion of meters to degrees? Need to think on that one

[dgboss]

Oh, I think I figured out one path forward, but we might want to dev chat next week. Right now when this warps, we don't specify the dimensions of the array. We can use height and width parameters in WarpOptions and set those to self.band.YSize and self.band.XSize. Also, we can keep the outputBounds and pass in outputBoundsSRS. Untested as of yet, but my brain is done.

[dgboss]

My brain is working again (sort of). Specifying the number of pixels for height and width might not be a good idea afterall. I think it could lead to non-square pixels when re-projecting.

[dgboss]

Closing to free up resources.

[sonarqubecloud]

Quality Gate passed

Issues
1 New issue
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
0.0% Duplication on New Code

See analysis details on SonarCloud