Pygeoprocessing error on warp_raster

Hi guys, I couldn't remember if you had a formal bug reporting method for Pygeoprocessing so I'l post here.

I have a problem when running the python version of pgp when it calls warp_raster. See the traceback below. I've tried several things including checking that it's finding my gdal data (it is) and giving it different option lists (doesn't help). Reverting to the previous version of warp_raster (i think from ~1.0) did work on the same inputs. Not time sensitive, but just fyi. Let me know if there's more info I should provide.



Error
Traceback (most recent call last):
  File "C:\Anaconda\envs\Projects\lib\unittest\case.py", line 59, in testPartExecutor
    yield
  File "C:\Anaconda\envs\Projects\lib\unittest\case.py", line 605, in run
    testMethod()
  File "C:\OneDrive\Projects\hazelbean\tests\test_spatial_utils.py", line 111, in test_clip_raster_by_vector
    hb.clip_raster_by_vector(self.global_1deg_raster_path, temp_path, self.two_polygon_shapefile_path, all_touched=True, ensure_fits=True)
  File "C:\OneDrive\Projects\hazelbean\hazelbean\spatial_utils.py", line 235, in clip_raster_by_vector
    gtiff_creation_options=gtiff_creation_options)
  File "C:\OneDrive\Projects\hazelbean\hazelbean\geoprocessing_extension.py", line 179, in align_and_resize_raster_stack_ensuring_fit
    gtiff_creation_options=option_list)
  File "C:\Anaconda\envs\Projects\lib\site-packages\pygeoprocessing\geoprocessing.py", line 1860, in warp_raster
    cutlineWhere=mask_vector_where_filter)
  File "C:\Anaconda\envs\Projects\lib\site-packages\osgeo\gdal.py", line 575, in Warp
    return wrapper_GDALWarpDestName(destNameOrDestDS, srcDSTab, opts, callback, callback_data)
  File "C:\Anaconda\envs\Projects\lib\site-packages\osgeo\gdal.py", line 3129, in wrapper_GDALWarpDestName
    return _gdal.wrapper_GDALWarpDestName(*args)
TypeError: in method 'wrapper_GDALWarpDestName', argument 4 of type 'GDALWarpAppOptions *'

Tagged:

Comments

  • RichRich Administrator, NatCap Staff
    Ooo lookit this! A pygeoprocessing issue!

    Can you post the line of code that's a call to `pygeoprocessing.warp_raster`? That looks like "File "C:\OneDrive\Projects\hazelbean\hazelbean\geoprocessing_extension.py", line 179"? 

    Sometimes when I see stuff like this it's as simple as there's a `unicode` string that needs to be cast to a `str` or a numpy numeric type that needs to be cast to a float... Now that I think about it, could you post as much as you think we'd need to recreate this error? I've been tolerating this privately, but maybe we should fix it internally to the `warp_raster` call.

  • RichRich Administrator, NatCap Staff
    And p.s. I'm happy to take this stuff on the forums, or you can also post an issue directly in the repo at https://bitbucket.org/natcap/pygeoprocessing/issues
  • Hey rich, I've attached my modified function that calls warp_raster and the testing function for calling it. While we're here, though, I thought I'd mention one more thing. This function was my half-assed way of fixing a problem with warp_raster whereby the resultant raster would shift upper-left to be alligned with the input, but this would cause the (sometimes undesired) feature of making the bottom right no longer be covered. My modification just adds an extra row and column in that case. Would there have been a more elegant way of doing this?

    base_raster_path_list = [self.global_ha_per_cell_5m]
    target_raster_path_list = [hb.temp('.tif', 'clip1', True)]
    resample_method_list = ['bilinear']
    target_pixel_size = get_raster_info(self.global_ha_per_cell_5m)['pixel_size']
    bounding_box_mode = 'intersection'
    base_vector_path_list = [self.two_poly_wgs84_aoi_path]
    raster_align_index = 0


    hb.align_and_resize_raster_stack_ensuring_fit(
    base_raster_path_list, target_raster_path_list, resample_method_list,
    target_pixel_size, bounding_box_mode, base_vector_path_list=base_vector_path_list,
    raster_align_index=raster_align_index, all_touched=True,
    gtiff_creation_options=hb.DEFAULT_GTIFF_CREATION_OPTIONS)





    def align_and_resize_raster_stack_ensuring_fit(
    base_raster_path_list, target_raster_path_list, resample_method_list,
    target_pixel_size, bounding_box_mode, base_vector_path_list=None,
    raster_align_index=None, ensure_fits=False, all_touched=False,
    gtiff_creation_options=hb.DEFAULT_GTIFF_CREATION_OPTIONS):
    """Generate rasters from a base such that they align geospatially.

    This function resizes base rasters that are in the same geospatial
    projection such that the result is an aligned stack of rasters that have
    the same cell size, dimensions, and bounding box. This is achieved by
    clipping or resizing the rasters to intersected, unioned, or equivocated
    bounding boxes of all the raster and vector input.

    Parameters:
    base_raster_path_list (list): a list of base raster paths that will
    be transformed and will be used to determine the target bounding
    box.
    target_raster_path_list (list): a list of raster paths that will be
    created to one-to-one map with `base_raster_path_list` as aligned
    versions of those original rasters.
    resample_method_list (list): a list of resampling methods which
    one to one map each path in `base_raster_path_list` during
    resizing. Each element must be one of
    "nearest|bilinear|cubic|cubic_spline|lanczos|mode".
    target_pixel_size (tuple): the target raster's x and y pixel size
    example: [30, -30].
    bounding_box_mode (string): one of "union", "intersection", or
    a list of floats of the form [minx, miny, maxx, maxy]. Depending
    on the value, output extents are defined as the union,
    intersection, or the explicit bounding box.
    base_vector_path_list (list): a list of base vector paths whose
    bounding boxes will be used to determine the final bounding box
    of the raster stack if mode is 'union' or 'intersection'. If mode
    is 'bb=[...]' then these vectors are not used in any calculation.
    raster_align_index (int): indicates the index of a
    raster in `base_raster_path_list` that the target rasters'
    bounding boxes pixels should align with. This feature allows
    rasters whose raster dimensions are the same, but bounding boxes
    slightly shifted less than a pixel size to align with a desired
    grid layout. If `None` then the bounding box of the target
    rasters is calculated as the precise intersection, union, or
    bounding box.
    gtiff_creation_options (list): list of strings that will be passed
    as GDAL "dataset" creation options to the GTIFF driver, or ignored
    if None.

    Returns:
    None
    """
    last_time = time.time()


    # make sure that the input lists are of the same length
    list_lengths = [
    len(base_raster_path_list), len(target_raster_path_list),
    len(resample_method_list)]
    if len(set(list_lengths)) != 1:
    raise ValueError(
    "base_raster_path_list, target_raster_path_list, and "
    "resample_method_list must be the same length "
    " current lengths are %s" % (str(list_lengths)))

    # we can accept 'union', 'intersection', or a 4 element list/tuple
    if bounding_box_mode not in ["union", "intersection"] and (
    not isinstance(bounding_box_mode, (list, tuple)) or
    len(bounding_box_mode) != 4):
    raise ValueError("Unknown bounding_box_mode %s" % (
    str(bounding_box_mode)))

    if ((raster_align_index is not None) and
    ((raster_align_index < 0) or
    (raster_align_index >= len(base_raster_path_list)))):
    raise ValueError(
    "Alignment index is out of bounds of the datasets index: %s"
    " n_elements %s" % (
    raster_align_index, len(base_raster_path_list)))

    raster_info_list = [
    get_raster_info(path) for path in base_raster_path_list]
    if base_vector_path_list is not None:
    vector_info_list = [
    get_vector_info(path) for path in base_vector_path_list]
    else:
    vector_info_list = []

    # get the literal or intersecting/unioned bounding box
    if isinstance(bounding_box_mode, (list, tuple)):
    target_bounding_box = bounding_box_mode
    else:
    # either intersection or union
    target_bounding_box = reduce(
    functools.partial(hb.merge_bounding_boxes, mode=bounding_box_mode),
    [info['bounding_box'] for info in
    (raster_info_list + vector_info_list)])

    if bounding_box_mode == "intersection" and (
    target_bounding_box[0] > target_bounding_box[2] or
    target_bounding_box[1] > target_bounding_box[3]):
    raise ValueError("The rasters' and vectors' intersection is empty "
    "(not all rasters and vectors touch each other).")

    if raster_align_index is not None:
    if raster_align_index >= 0:
    # bounding box needs alignment
    align_bounding_box = (
    raster_info_list[raster_align_index]['bounding_box'])
    align_pixel_size = (
    raster_info_list[raster_align_index]['pixel_size'])
    # adjust bounding box so lower left corner aligns with a pixel in
    # raster[raster_align_index]
    target_rc = [
    math.ceil((target_bounding_box[2] - target_bounding_box[0]) / target_pixel_size[0]),
    math.floor((target_bounding_box[3] - target_bounding_box[1]) / target_pixel_size[1])
    ]

    original_bounding_box = list(target_bounding_box)

    for index in [0, 1]:
    n_pixels = int((target_bounding_box[index] - align_bounding_box[index]) / float(align_pixel_size[index]))

    target_bounding_box[index] = align_bounding_box[index] + (n_pixels * align_pixel_size[index])
    target_bounding_box[index + 2] = target_bounding_box[index] + target_rc[index] * target_pixel_size[index]

    if ensure_fits:
    # This addition to the core geoprocessing code was to fix the case where the alignment moved the target tif
    # up and to the left, but in a way that then trunkated 1 row/col on the bottom right, causing wrong-shape
    # raster_math errors.z
    if original_bounding_box[2] > target_bounding_box[2]:
    target_bounding_box[2] += target_pixel_size[0]

    if original_bounding_box[3] > target_bounding_box[3]:
    target_bounding_box[3] -= target_pixel_size[1]


    option_list = list(gtiff_creation_options)
    if all_touched:
    option_list.append("ALL_TOUCHED=TRUE")

    for index, (base_path, target_path, resample_method) in enumerate(zip(
    base_raster_path_list, target_raster_path_list,
    resample_method_list)):
    last_time = gp._invoke_timed_callback(
    last_time, lambda: L.info(
    "align_dataset_list aligning dataset %d of %d",
    index, len(base_raster_path_list)), hb.LOGGING_PERIOD)
    option_list = []

    ## My replacement call to the older version
    # hb.warp_raster_HAZELBEAN_REPLACEMENT(
    # base_path, target_pixel_size,
    # target_path, resample_method,
    # target_bb=target_bounding_box,
    # gtiff_creation_options=option_list)

    # # PGP Replacement: warp_was give gdal error.
    gp.warp_raster(
    base_path, target_pixel_size,
    target_path, resample_method,
    target_bb=target_bounding_box,
    gtiff_creation_options=option_list)
  • RichRich Administrator, NatCap Staff
    I think I get what you're trying to do... Let me see. It's true that in warp raster we don't extend a fractional pixel into a whole pixel. We used to do that, and we'd get into even crazier cases where there'd be baaarely an overlap and we'd have ourselves a new row and column in a raster.

    So the whole fix to this in PGP is to be consistent that we will not create fractional pixels. In turn, to make sure everything is all aligned, the intent of `align_and_resize_raster_stack` is send the ENTIRE raster stack all at once. That way you don't have to try to hope that whatever algorithm we chose will kill your manicured raster stack. 

    BUT we also get that not everyone would want to do this and in some cases they already know their raster stack is aligned but they want to just resize/warp one raster. So that's why we added a `target_bb` optional argument in `warp_raster` that can be gotten from a call to `get_raster_info(path)['bounding_box']`.  

    From your first code snippet, I can't quite tell what you're trying to do. It looks like you're passing one raster, but you're also using the raster_align_index. If you actually have a whole raster list to process, would you consider aligning the whole thing as a stack? 

    BUT ALSO, I'm gonna guess this is this raster you have that you pre-calculated the areas of cells on a WGS84 projection, right? If it's helpful, in the most recent version of PyGeoprocessing we added a feature to `raster_calculator` to allow you to pass 1 and 2 dimensional arrays that would otherwise project to the shape of the input rasters. I have used this to include a dynamic "area in pixel" calculation by first calculating the area of a pixel in a 1 dimensional column for all the latitudes that are relevant to my area. From there I pass it directly into raster calculator and inside the raster calculator op it gets treated as any other 2d input array. If you're  curious, check out `calc_pop_count` here: https://bitbucket.org/natcap/ipbes-pollination-analysis/src/da484864278dc7a95496bfac54bbe734e46fa25c/ipbes_pollination.py#lines-1754


    But that's all to say I can't see what the issue might be for your original issue. Could you maybe try to make sure your `self.global_ha_per_cell_5m`, `hb.temp('.tif', 'clip1', True)` and `self.two_poly_wgs84_aoi_path` are cast to `str` before you make the call? And if that doesn't work, since you have your own working copy of pygeoprocessing, can you dump every argument to the call

    gp.warp_raster(
                base_path, target_pixel_size,
                target_path, resample_method,
                target_bb=target_bounding_box,
                gtiff_creation_options=option_list)

    Or enter the pdb debugger there and see if any parameters look unusual?
Sign In or Register to comment.