Raster-based Spatial Analysis for Python

Addresses #125, #124

Using the circle kernel to extend to an annulus kernel. Implementation excludes the cells on the bound of the inner radius of the annulus.

Also implements the ability to return the raw z-scores from hotspots.

versioneer is deprecated, does not work properly with pep517/518, and the source distribution tarball on PyPI is broken without the versioneer.py file in it. See https://github.com/conda-forge/staged-recipes/pull/12177

This PR implements setuptools_scm instead of versioneer, modernize the Travis-CI a little bit and implement a tarball check to avoid regressions, and fix the tarball.

Ping @rsignell-usgs who requested the conda-forge package.

Needed in https://github.com/conda-forge/staged-recipes/pull/12177

The new code uses less context switching and data transfers, minimizes communication overhead and removes the cupy dependency. I've made it run by default, if cuda is available. Otherwise if falls back to numpy

This change addresses issue #541

Given a 2D raster image, a set of target pixels in the input raster, the proximity function computes a raster of proximity to indicate the distance from each pixel to the nearest target pixel. The distance metric to be used can be one of the following: Euclidean, Great-Circle, or Manhattan. In the most naive way, we can calculate distances from each pixel to all target pixels and find the closest target. This would take m*n*t calculations, where m*n is the size of the raster, and t is the number of targets.

The current implementation of xrspatial.proximity function is ported from GDAL with some modifications to make it work with xarray.DataArray. This notebook shows how to use the function. To keep it simple, let's consider the problem at a 2D array level instead of a 2D xarray DataArray. The algorithm can be described as follows:

Inputs:

• Raster image I, a height x width 2D array
• Set of target pixels T. All target pixels are in I
• Distance metric d()

Output: Proximity raster P where P[i, j] is the distance from cell (i, j) to its nearest target pixel

Idea: Use dynamic programming to identify the nearest target pixel of a pixel based on the nearest target pixels of its 3x3 neighborhood window.

Detail implementation:

• Let Nx be a 1d array of width elements: Nx[j] is the x-position in pixel space of nearest target pixel of pixel (i, j). Nx[j] can have a value in [0, width-1]
• Let Ny be a 1d array of width elements: Ny[j] is the y-position in pixel space of nearest target pixel of pixel (i, j). Ny[j] can have a value in [0, height-1]. Values of Nx and Ny will be updated row by row.

1. Initially,

• set Nx[j] = -1 for all j.
• set Ny[j] = -1 for all j.
• for all (i, j), set P[i, j] = 0 if cell (i, j) is a target pixel, P[i, j] = infinity otherwise.
• note that distance d( (i1, j1), (i2, j2) ) = infinity if any of the inputs i1, j1, i2, j2 equals to -1 (i.e, invalid pixel)

2. Traverse the image row by row from top to bottom:

• traverse each row from left to right
• traverse each row from right to left

3. Reset Nx[j] = -1 and Ny[j] = -1 for all j.
4. Traverse the image in reverse order from bottom to top:

• traverse each row from right to left
• traverse each row from left to right

The formula to update P[i, j] and Ny[j] and Nx[j] at each cell (i, j)

P[i, j] = min(
    P[i, j],                   
    d( (Ny[j], Nx[j]), (i, j) ),     # Are we nearer to the closest target to the above (below) pixel?
    d( (Ny[j-1], Nx[j-1]), (i, j) ), # Are we nearer to the closest target to the left (right) pixel?
    d( (Ny[j+1], Nx[j+1]), (i, j) ), # Are we nearer to the closest target to the top right (bottom left) pixel?
)

Update Ny[j] and Nx[j] accordingly:

• Ny[j], Nx[j] = i, j if P[i, j] = 0
• Ny[j], Nx[j] = Ny[j-1], Nx[j-1] if P[i, j] is updated as d( (Ny[j-1], Nx[j-1]), (i, j) )
• Ny[j], Nx[j] = Ny[j+1], Nx[j+1] if P[i, j] is updated as d( (Ny[j+1], Nx[j+1]), (i, j) )

We're looking for an implementation for the distributed version of the Proximity function that works with Dask. Explicit questions are listed as:

1. Currently, the calculations are performed sequentially, how to parallelize them?
2. Dividing our 2D input data into a set of smaller chunks, how to compute proximity distance grid chunk by chunk? The primary question is, how to determine the nearest target pixel of all pixels in a chunk when the target pixels can be outside the chunk?

Describe the bug Running zonal_stats and passing in the argument nodata_values = -99999 results in a ValueError: All arrays must be of the same length stack trace with method and line number zonal.py line 519 zonal.py 335 python3.7/site-packages/pandas/core/frame.py 614 python3.7/site-packages/pandas/core/internals/construction.py 465

When i run zonal_stats without the nodata_values argument it runs fine and computes

nodata_values = -99999 or nodata_values = -99999.0

Expected behavior Would expect the method zonal_stats to return statistics computed ignoring pixels with nodata_values

Screenshots Unable to provide screenshots, running this in a classified environment

Desktop (please complete the following information): Running this in Databricks 7.3LTS with Python 3.7.5, xarray-spatial 0.3.0

Describe the bug Cannot install the latest version of xarray-spatial with an up-to-date stack (numpy>1.20).

Expected behavior Install the latest version of xarray-spatial through conda

Screenshots I'm trying to update the gds_env and trying an install that pins xarray-spatial to its latest release returns the following error:

package xarray-spatial-0.2.9-pyhd8ed1ab_0 requires numpy >=1.7,<=1.20, but none of the providers can be installed

The numpy version it is installing (for compatibility with other libraries) is 1.21.2. Is there any reason why it pins to 1.20 as the most recent version?

Added reprojection with one blank function to test optional modules install. Can be extended to other modules as well. Install with pip install -e . --install-option='--reprojection' to get reprojection module.

Changed versioning scheme. The scm versioning in setup py was not working correctly. I substituted from param (as in datashader) to generate versions. Conda versions are fixed at 0.1.6 right now.

Now that numba 0.55 has been released with support for numpy 1.21 and python 3.10, we no longer need to pin to a max numpy version ourselves.

Fixes issue #565.

This PR brings all the terrain generation related development and optimization. In brief:

• adds a fast_perlin.py module that contains a fast_perlin() method that can be used to generate perlin noise. It offers higher performance compared to the previous implementation of perlin noise and has a Numpy, CuPy and Dask implementation. Note that the interface is now slightly different. The user needs to pass an xr.DataArray (can be based on numpy, cupy or dask array) instead of width, height.
• adds a fast_terrain.py module that contains a generate_fast_terrain() method that can be used to generate a randomized terrain. It is faster compare to the existing implementation of terrain generation. Uses the fast_perlin() method and has a Numpy, CuPy and Dask implementation. Note that the interface is now slightly different. The user needs to pass an xr.DataArray (can be based on numpy, cupy or dask array) instead of a Canvas.
• Adds two jupyter notebooks in the examples directory, one for each of the two newly added methods.

Hey the last few releases need the CHANGELOG.md information copied over. to the github Release for each version (e.g. 0.1.4, 0.1.5, 0.1.6). Check those out, the content should be fine the way it is.

I'm testing xarray-spatial against QGIS and seeing that when running hillshade on the same input data, the results by xarray-spatial and GDAL/QGIS are very different. Source code for QGIS hillshade can be found at: https://github.com/qgis/QGIS/blob/master/python/plugins/processing/algs/gdal/hillshade.py. This needs more research to see how the algorithm was implemented in QGIS to clearly understand the difference between the 2 libraries.

Input data:

array([[      nan,       nan,       nan,       nan,       nan,       nan],
       [704.237  , 242.24084, 429.3324 , 779.8816 , 193.29506, 984.6926 ],
       [226.56795, 815.7483 , 290.6041 ,  76.49687, 820.89716,  32.27882],
       [344.8238 , 256.34998, 806.8326 , 602.0442 , 721.1633 , 496.95636],
       [185.43515, 834.10425, 387.0871 , 716.0262 ,  49.61273, 752.95483],
       [302.4271 , 151.49211, 442.32797, 358.4702 , 659.8187 , 447.1241 ],
       [148.04834, 819.2133 , 468.97913, 977.11694, 597.69666, 999.14185],
       [268.1575 , 625.96466, 840.26483, 448.28333, 859.2699 , 528.04095]],
      dtype=float32)

xarray-spatial hillshade:

array([[       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan, 0.75030494, 0.06941041, 0.90643436, 0.15474272,            nan],
       [       nan, 0.80836594, 0.72366774, 0.14052185, 0.774778  ,            nan],
       [       nan, 0.93396175, 0.7071851 , 0.42872226, 0.9455124 ,            nan],
       [       nan, 0.85551083, 0.6819584 , 0.46013114, 0.23561102,            nan],
       [       nan, 0.41484872, 0.3213355 , 0.5821109 , 0.21879822,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan]],
 dtype=float32)

QGIS/GDAL hillshade

array([[  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [107.84468 ,  70.09885 ,   0.      ,  17.661407,   0.      ,          0.      ],
       [ 80.06987 ,  71.644684,   0.      ,   0.      ,   0.      ,          0.      ],
       [ 85.574615, 106.36669 ,  96.23605 ,  28.27108 ,  90.29079 ,         85.07072 ],
       [ 81.44522 ,  77.092354,   8.479876,   0.      ,   0.      ,          0.      ],
       [ 62.541145,   2.647696,   0.      ,   0.      ,   0.      ,          6.515689],
       [ 74.07955 ,  78.71434 ,   0.      ,  84.590744,  34.814816,         44.81609 ]],
 dtype=float32)

Many GIS data is collected frequently. For example, Sentinel-2 is captured every 5 days. Having tools that support studying data over time can help in detecting the changes, as well as identifying the trends to predict future. Time series analysis will be our next step. Let's do some research about how to do time series hotspots analysis.

https://pro.arcgis.com/en/pro-app/latest/tool-reference/space-time-pattern-mining/emerginghotspots.htm

https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-statistics/space-time-analysis.htm

One of the 2023 focus is to strengthen our GPU support to xarray-spatial tools. Let's review again and address which needs improvement and what we can do.

To evaluate the reliability of xarray-spatial, let's test the package and compare the results against QGIS/GDAL.

• [x] aspect
• [x] slope
• [x] multispectral tools
• [x] proximity
• [x] zonal stats

Fixes #22, #104, #182