Numba/CUDA projection kernels for reproject, README update

README.md

@@ -79,13 +79,15 @@
 
 :fast_forward: Scalable with [Dask](http://dask.pydata.org)
 
+:desktop_computer: GPU-accelerated with [CuPy](https://cupy.dev/) and [Numba CUDA](https://numba.readthedocs.io/en/stable/cuda/index.html)
+
 :confetti_ball: Free of GDAL / GEOS Dependencies
 
 :earth_africa: General-Purpose Spatial Processing, Geared Towards GIS Professionals
 
 -------
 
-Xarray-Spatial implements common raster analysis functions using Numba and provides an easy-to-install, easy-to-extend codebase for raster analysis.
+Xarray-Spatial is a Python library for raster analysis built on xarray. It has 100+ functions for surface analysis, hydrology (D8, D-infinity, MFD), fire behavior, flood modeling, multispectral indices, proximity, classification, pathfinding, and interpolation. Functions dispatch automatically across four backends (NumPy, Dask, CuPy, Dask+CuPy). A built-in GeoTIFF/COG reader and writer handles raster I/O without GDAL.
 
 ### Installation
 ```bash
@@ -119,9 +121,9 @@ In all the above, the command will download and store the files into your curren
 
 `xarray-spatial` grew out of the [Datashader project](https://datashader.org/), which provides fast rasterization of vector data (points, lines, polygons, meshes, and rasters) for use with xarray-spatial.
 
-`xarray-spatial` does not depend on GDAL / GEOS, which makes it fully extensible in Python but does limit the breadth of operations that can be covered.  xarray-spatial is meant to include the core raster-analysis functions needed for GIS developers / analysts, implemented independently of the non-Python geo stack.
+`xarray-spatial` does not depend on GDAL or GEOS. Raster I/O, reprojection, compression codecs, and coordinate handling are all pure Python and Numba -- no C/C++ bindings anywhere in the stack.
 
-Our documentation is still under construction, but [docs can be found here](https://xarray-spatial.readthedocs.io/en/latest/).
+[API reference docs](https://xarray-spatial.readthedocs.io/en/latest/) and [33+ user guide notebooks](examples/user_guide/) cover every module.
 
 #### Raster-huh?
 
@@ -210,9 +212,46 @@ write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
 
 | Name | Description | Source | NumPy xr.DataArray | Dask xr.DataArray | CuPy GPU xr.DataArray | Dask GPU xr.DataArray |
 |:----------:|:------------|:------:|:----------------------:|:--------------------:|:-------------------:|:------:|
-| [Reproject](xrspatial/reproject/__init__.py) | Reprojects a raster to a new CRS using an approximate transform and numba JIT resampling | Standard (inverse mapping) | ✅️ | ✅️ | ✅️ | ✅️ |
+| [Reproject](xrspatial/reproject/__init__.py) | Reprojects a raster to a new CRS with Numba JIT / CUDA coordinate transforms and resampling | Standard (inverse mapping) | ✅️ | ✅️ | ✅️ | ✅️ |
 | [Merge](xrspatial/reproject/__init__.py) | Merges multiple rasters into a single mosaic with configurable overlap strategy | Standard (mosaic) | ✅️ | ✅️ | 🔄 | 🔄 |
 
+Built-in Numba JIT and CUDA projection kernels bypass pyproj for common CRS pairs. Other CRS pairs fall back to pyproj automatically.
+
+| Projection | EPSG examples | CPU Numba | CUDA GPU |
+|:-----------|:-------------|:---------:|:--------:|
+| Web Mercator | 3857 | ✅️ | ✅️ |
+| UTM / Transverse Mercator | 326xx, 327xx, State Plane | ✅️ | ✅️ |
+| Ellipsoidal Mercator | 3395 | ✅️ | ✅️ |
+| Lambert Conformal Conic | 2154, State Plane | ✅️ | ✅️ |
+| Albers Equal Area | 5070 | ✅️ | ✅️ |
+| Cylindrical Equal Area | 6933 | ✅️ | ✅️ |
+| Sinusoidal | MODIS grids | ✅️ | ✅️ |
+| Lambert Azimuthal Equal Area | 3035, 6931, 6932 | ✅️ | ✅️ |
+| Polar Stereographic | 3031, 3413, 3996 | ✅️ | ✅️ |
+
+**Reproject performance** (end-to-end, bilinear, vs rioxarray):
+
+| Transform | 1024x1024 | | 4096x4096 | |
+|:---|---:|---:|---:|---:|
+| | xrspatial | rioxarray | xrspatial | rioxarray |
+| WGS84 -> UTM 33N | 33ms | 72ms (2.2x) | 627ms | 1.09s (1.7x) |
+| WGS84 -> Web Mercator | 16ms | 44ms (2.9x) | 526ms | 741ms (1.4x) |
+| WGS84 -> Albers CONUS | 72ms | 196ms (2.7x) | 649ms | 1.78s (2.7x) |
+| WGS84 -> LAEA Europe | 47ms | 74ms (1.6x) | 677ms | 1.03s (1.5x) |
+| WGS84 -> Polar Stere S | 34ms | 580ms (17x) | 839ms | 9.13s (11x) |
+
+Times include coordinate transform + bilinear resampling. Speedup in parentheses is rioxarray/xrspatial. The Polar Stereographic advantage comes from rioxarray computing a much larger output grid for the same input extent.
+
+**Merge performance** (4 overlapping same-CRS tiles, vs rioxarray):
+
+| Tile size | xrspatial | rioxarray | Speedup |
+|:---|---:|---:|---:|
+| 512x512 | 11ms | 50ms | **4.5x** |
+| 1024x1024 | 82ms | 125ms | **1.5x** |
+| 2048x2048 | 347ms | 604ms | **1.7x** |
+
+Same-CRS tiles skip reprojection entirely and are placed by direct coordinate alignment.
+
 -------
 
 ### **Utilities**
@@ -460,29 +499,29 @@ write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
 Importing `xrspatial` registers an `.xrs` accessor on DataArrays and Datasets, giving you tab-completable access to every spatial operation:
 
 ```python
-import xrspatial
-from xrspatial.geotiff import read_geotiff
+import xrspatial as xrs
+from xrspatial.geotiff import read_geotiff, write_geotiff
 
 # Read a GeoTIFF (no GDAL required)
 elevation = read_geotiff('dem.tif')
 
-# Surface analysis — call operations directly on the DataArray
+# Surface analysis
 slope = elevation.xrs.slope()
 hillshaded = elevation.xrs.hillshade(azimuth=315, angle_altitude=45)
 aspect = elevation.xrs.aspect()
 
+# Reproject and write as a Cloud Optimized GeoTIFF
+dem_wgs84 = elevation.xrs.reproject(target_crs='EPSG:4326')
+write_geotiff(dem_wgs84, 'output.tif', cog=True)
+
 # Classification
 classes = elevation.xrs.equal_interval(k=5)
 breaks = elevation.xrs.natural_breaks(k=10)
 
 # Proximity
 distance = elevation.xrs.proximity(target_values=[1])
 
-# Multispectral — call on the NIR band, pass other bands as arguments
-nir = xr.DataArray(np.random.rand(100, 100), dims=['y', 'x'])
-red = xr.DataArray(np.random.rand(100, 100), dims=['y', 'x'])
-blue = xr.DataArray(np.random.rand(100, 100), dims=['y', 'x'])
-
+# Multispectral
 vegetation = nir.xrs.ndvi(red)
 enhanced_vi = nir.xrs.evi(red, blue)
 ```
@@ -503,14 +542,14 @@ ndvi_result = ds.xrs.ndvi(nir='band_5', red='band_4')
 
 ##### Function Import Style
 
-All operations are also available as standalone functions if you prefer explicit imports:
+All operations are also available as standalone functions:
 
 ```python
-from xrspatial import hillshade, slope, ndvi
+import xrspatial as xrs
 
-hillshaded = hillshade(elevation)
-slope_result = slope(elevation)
-vegetation = ndvi(nir, red)
+hillshaded = xrs.hillshade(elevation)
+slope_result = xrs.slope(elevation)
+vegetation = xrs.ndvi(nir, red)
 ```
 
 Check out the user guide [here](/examples/user_guide/).