This forum is shutting down! Please post new discussions at community.naturalcapitalproject.org

Pygeoprocessing Vectoring Operation

Hello,

This is in reference to PyGeoprocessing (https://bitbucket.org/richpsharp/pygeoprocessing). I was wondering if you could give a simple example of a vectorized and non-vectorized operation for pygeoprocessing.geoprocessing.vecotrize_datasets. For example, increasing a raster by a constant while preserving nodata values.

Thank you,
Martin

Comments

  • RichRich Administrator, NatCap Staff
    Hi Martin, see below for a complete example of using a non-vectorized operation for vectorize_datasets.  I won't show the vectorized form since it only exists for backwards compatibility on a couple remaining InVEST models.  That option will be removed entirely in the future and is already deprecated.  Though if you were dying for examples, you could look at the source code in natcap.invest.scenario_generator :)

    """Complete, but trival example of using vectorize_datasets to increment values
    in a raster by a constant"""
    import sys

    import gdal
    import numpy
    import pygeoprocessing

    def add_one(input_raster_filename, constant, output_raster_filename):
        """Adds `constant` to the values in the raster found at
        `input_raster_filename` and places the result in `output_raster_filename`
        which will be a 32 bit floating point raster of the same projection and
        shape as the raster located at `input_raster_filename.  If
        `output_raster_filename` exists, it is overwritten with the result. The
        nodata values in input_raster_filename are ignored when incrementing the
        value."""

        out_pixel_size = pygeoprocessing.geoprocessing.get_cell_size_from_uri(
            input_raster_filename)
        nodata = pygeoprocessing.geoprocessing.get_nodata_from_uri(
            input_raster_filename)

        def _add_one_op(in_array):
            """increments input array by constant unless it is a nodata value"""
            return numpy.where(in_array != nodata, in_array + constant, nodata)

        pygeoprocessing.geoprocessing.vectorize_datasets(
            [input_raster_filename], _add_one_op, output_raster_filename,
            gdal.GDT_Float32, nodata, out_pixel_size, "intersection",
            vectorize_op=False)

    if __name__ == '__main__':
        add_one(sys.argv[1], float(sys.argv[2]), sys.argv[3])

    If the code above were in a file called `vectorize_datasets.py` you could invoke it at the command line like this:

    python vectorize_datasets.py edge_carbon_lu_sample_map.tif 123.5 out.tif

    Where `edge_carbon_lu_sample_map.tif`  would be a GIS file that already exists, and `123.5` is your constant and `out.tif` is your output file.
  • Thanks Rich!
Sign In or Register to comment.