Added warp, ranges.

.github/workflows/CI.yml

@@ -38,7 +38,13 @@ jobs:
             ${{ runner.os }}-test-
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
+      - name: Set PROJ_CURL_CA_BUNDLE environment variable
+        if: runner.os == 'Linux'
+        run: |
+          echo "PROJ_CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt" >> $GITHUB_ENV
       - uses: julia-actions/julia-runtest@v1
+        env: 
+          PROJ_NETWORK: "ON"
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v2
         with:

CHANGELOG.md

@@ -7,6 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+## [0.8.0] - 2023-03-27
+- Added `warp` for warping GeoArrays.
+- Added `ranges` for returning the x and y `StepRange` of coordinates.
+- Replaced `equals` with `Base.isequal` and made sure to compare the AffineMap only approximately to account for floating point precision.
+- `coords(ga)` now returns an iterator. Apply `collect` on it for the old behaviour.
+- `indices` now returns a `CartesianIndex` instead of `i, j`. Call `.I` on it for the old behaviour.
+
 ## [0.7.13] - 2023-01-12
 - Added convert, affine!
 - Added broadcast for GeoArray, so `ga .+ 1` isa `GeoArray`
@@ -39,7 +46,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Fixed iterate specification so `sum` on a GeoArray is correct
 
-[unreleased]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.13...HEAD
+[unreleased]: https://github.com/evetion/GeoArrays.jl/compare/v0.8.0...HEAD
+[0.8.0]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.13...v0.8.0
 [0.7.13]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.12...v0.7.13
 [0.7.12]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.11...v0.7.12
 [0.7.11]: https://github.com/evetion/GeoArrays.jl/compare/v0.7.10...v0.7.11

Project.toml

@@ -1,11 +1,13 @@
 name = "GeoArrays"
 uuid = "2fb1d81b-e6a0-5fc5-82e6-8e06903437ab"
 authors = ["Maarten Pronk <[email protected]>"]
-version = "0.7.13"
+version = "0.8.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoStatsBase = "323cb8eb-fbf6-51c0-afd0-f8fba70507b2"
@@ -16,6 +18,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 [compat]
 ArchGDAL = "0.7 - 0.9, 0.10"
 CoordinateTransformations = "0.5 - 0.6"
+DataAPI = "1"
+Extents = "0.1"
 GeoFormatTypes = "0.4"
 GeoInterface = "1"
 GeoStatsBase = "0.21 - 0.30"

README.md

@@ -11,7 +11,7 @@ This packages takes its inspiration from Python's [rasterio](https://github.com/
 ## Installation
 
 ```julia
-(v1.7) pkg> add GeoArrays
+(v1.8) pkg> add GeoArrays
 ```
 
 ## Examples
@@ -85,30 +85,35 @@ In case there is missing data, the type will be a `Union{Missing, T}`. To conver
 
 ```julia
 # Find coordinates by index
-julia> GeoArrays.coords(geoarray, [1,1])
+julia> GeoArrays.coords(geoarray, (1,1))
 2-element StaticArrays.SArray{Tuple{2},Float64,1,2}:
  440720.0
       3.75132e6
 ```
 
-All coordinates (tuples) are obtained when omitting the index parameter.
+All coordinates (tuples) are obtained as generator when omitting the index parameter.
 
 ```julia
 # Find all coordinates
-julia> GeoArrays.coords(geoarray)
-101×101 Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}:
+julia> collect(GeoArrays.coords(geoarray))
+101×101 Matrix{StaticArraysCore.SVector{2, Float64}}:
  [440720.0, 3.75132e6]  [440720.0, 3.75126e6]  [440720.0, 3.7512e6] ...
  ...
 ```
 
+Similarly, one can find the coordinates ranges of a GeoArray
+
+```julia
+julia> x, y = GeoArrays.ranges(geoarray)
+(440750.0:60.0:446690.0, 3.75129e6:-60.0:3.74535e6)
+```
+
 The operation can be reversed, i.e. row and column index can be computed from coordinates with the `indices` function.
 
 ```julia
 # Find index by coordinates
 julia> indices(geoarray, [440720.0, 3.75132e6])
-2-element StaticArrays.SArray{Tuple{2},Int64,1,2}:
- 1
- 1
+CartesianIndex(1, 1)
 ```
 
 ### Manipulation
@@ -127,6 +132,19 @@ julia> GeoArray(rand(5,5,1)) - GeoArray(rand(5,5,1))
 5x5x1 Array{Float64,3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
 ```
 
+One can also warp an array, using GDAL behind the scenes.
+For example, we can vertically transform from the ellipsoid
+to the EGM2008 geoid using EPSG code 3855.
+Note that the underlying PROJ library needs to find the geoidgrids,
+so if they're not available locally, one needs to set `ENV["PROJ_NETWORK"] = "ON"`
+as early as possible, ideally before loading GeoArrays.
+```julia
+ga = GeoArray(zeros((360, 180)))
+bbox!(ga, (min_x=-180, min_y=-90, max_x=180, max_y=90))
+crs!(ga, GeoFormatTypes.EPSG(4979))  # WGS83 in 3D (reference to ellipsoid)
+ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
+```
+
 ### Nodata filling
 GeoArrays with missing data can be filled with the [`fill!`](@ref) function.
 
@@ -194,3 +212,8 @@ julia> plot(ga_sub)
 
 ### Profile
 You can sample the values along a line in a GeoArray with `profile(ga, linestring)`. The linestring can be any geometry that supports [GeoInterface.jl](https://github.com/JuliaGeo/GeoInterface.jl/).
+
+
+## Alternatives
+GeoArrays.jl was written to quickly save a geospatial Array to disk. Its functionality mimics `rasterio` in Python. If one requires more features---such as rasterization or zonal stats---which also work on NetCDF files, [Rasters.jl](https://github.com/rafaqz/Rasters.jl/) is a good alternative. Its functionality is more like `(rio)xarray` in Python.
+

docs/src/index.md

@@ -1,20 +1,20 @@
 ```@meta
 CurrentModule = GeoArrays
 ```
-[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://evetion.github.io/GeoArrays.jl/dev) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://evetion.github.io/GeoArrays.jl/stable) [![CI](https://github.com/evetion/GeoArrays.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/evetion/GeoArrays.jl/actions/workflows/CI.yml) [![codecov](https://codecov.io/gh/evetion/GeoArrays.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/evetion/GeoArrays.jl)
+[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://evetion.github.io/GeoArrays.jl/dev) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://evetion.github.io/GeoArrays.jl/stable) [![CI](https://github.com/evetion/GeoArrays.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/evetion/GeoArrays.jl/actions/workflows/CI.yml) [![codecov](https://codecov.io/gh/evetion/GeoArrays.jl/branch/master/graph/badge.svg?token=pXywTRCftQ)](https://codecov.io/gh/evetion/GeoArrays.jl)
 
 # GeoArrays
 
 Simple geographical raster interaction built on top of [ArchGDAL](https://github.com/yeesian/ArchGDAL.jl/), [GDAL](https://github.com/JuliaGeo/GDAL.jl) and [CoordinateTransformations](https://github.com/FugroRoames/CoordinateTransformations.jl).
 
-A [`GeoArray`](@ref) is an AbstractArray, an AffineMap for calculating coordinates based on the axes and a CRS definition to interpret these coordinates into in the real world. It's three dimensional and can be seen as a stack (3D) of 2D geospatial rasters (bands), the dimensions are :x, :y, and :bands. The AffineMap and CRS (coordinates) only operate on the :x and :y dimensions.
+A GeoArray is an AbstractArray, an AffineMap for calculating coordinates based on the axes and a CRS definition to interpret these coordinates into in the real world. It's three dimensional and can be seen as a stack (3D) of 2D geospatial rasters (bands), the dimensions are :x, :y, and :bands. The AffineMap and CRS (coordinates) only operate on the :x and :y dimensions.
 
 This packages takes its inspiration from Python's [rasterio](https://github.com/mapbox/rasterio).
 
 ## Installation
 
 ```julia
-(v1.7) pkg> add GeoArrays
+(v1.8) pkg> add GeoArrays
 ```
 
 ## Examples
@@ -27,7 +27,7 @@ Load the `GeoArrays` package.
 julia> using GeoArrays
 ```
 
-Read a GeoTIFF file, using [`GeoArrays.read`](@ref) and display its information, i.e. AffineMap and projection (CRS).
+Read a GeoTIFF file and display its information, i.e. AffineMap and projection (CRS).
 
 ```julia
 # Read TIF file
@@ -45,19 +45,19 @@ GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), "PROJCS[\
 ```
 
 ### Writing to GeoTIFF
-
-Create a random `GeoArray` and [`GeoArrays.write`](@ref) it to a GeoTIFF file. Setting the bounding box and the crs can be set by using [`bbox!`](@ref) and [`epsg!`](@ref).
+Create a random `GeoArray` and write it to a GeoTIFF file.
 
 ```julia
 # Create, reference and write a TIFF
 julia> ga = GeoArray(rand(100,200))
 julia> bbox!(ga, (min_x=2., min_y=51., max_x=5., max_y=54.))  # roughly the Netherlands
 julia> epsg!(ga, 4326)  # in WGS84
 julia> GeoArrays.write("test.tif", ga)
+# Or write it with compression and tiling
+julia> GeoArrays.write("test_compressed.tif", ga; options=Dict("TILED"=>"YES", "COMPRESS"=>"ZSTD"))
 ```
 
 ### Streaming support
-
 The package supports streaming reading.
 
 ```julia
@@ -68,10 +68,9 @@ julia> @time ga = GeoArrays.read(fn, masked=false)
 ```
 
 ### Reading bands
-
 GeoTIFFs can be large, with several bands, one can read.
 
-When working with large rasters, e.g. with satellite images that can be GB in size, it is useful to be able to read only one band (or a selection of them) to `GeoArray`. When using [`read`](@ref), one can specify the band.
+When working with large rasters, e.g. with satellite images that can be GB in size, it is useful to be able to read only one band (or a selection of them) to `GeoArray`. When using `read`, one can specify the band.
 
 ```julia
 # Get file
@@ -85,38 +84,42 @@ julia> ga_band = GeoArrays.read(fn, masked=false, band=2)
 In case there is missing data, the type will be a `Union{Missing, T}`. To convert to a GeoArray without `missing`, you can call `coalesce(ga, value_to_replace_missing)`.
 
 ### Using coordinates
-
-`GeoArray` has geographical coordinates for all array elements (pixels). They can be retrieved with the [`coords`](@ref) function.
+`GeoArray`s have geographical coordinates for all array elements (pixels). They can be retrieved with the `GeoArrays.coords` function.
 
 ```julia
 # Find coordinates by index
-julia> coords(geoarray, [1,1])
+julia> GeoArrays.coords(geoarray, (1,1))
 2-element StaticArrays.SArray{Tuple{2},Float64,1,2}:
  440720.0
       3.75132e6
 ```
 
-All coordinates (tuples) are obtained when omitting the index parameter.
+All coordinates (tuples) are obtained as generator when omitting the index parameter.
 
 ```julia
 # Find all coordinates
-julia> coords(geoarray)
-101×101 Array{StaticArrays.SArray{Tuple{2},Float64,1,2},2}:
+julia> collect(GeoArrays.coords(geoarray))
+101×101 Matrix{StaticArraysCore.SVector{2, Float64}}:
  [440720.0, 3.75132e6]  [440720.0, 3.75126e6]  [440720.0, 3.7512e6] ...
  ...
 ```
 
-The operation can be reversed, i.e. row and column index can be computed from coordinates with the [`indices`](@ref) function.
+Similarly, one can find the coordinates ranges of a GeoArray
+
+```julia
+julia> x, y = GeoArrays.ranges(geoarray)
+(440750.0:60.0:446690.0, 3.75129e6:-60.0:3.74535e6)
+```
+
+The operation can be reversed, i.e. row and column index can be computed from coordinates with the `indices` function.
 
 ```julia
 # Find index by coordinates
 julia> indices(geoarray, [440720.0, 3.75132e6])
-2-element StaticArrays.SArray{Tuple{2},Int64,1,2}:
- 1
- 1
+CartesianIndex(1, 1)
 ```
-### Manipulation
 
+### Manipulation
 Basic `GeoArray` manipulation is implemented, e.g. translation.
 ```julia
 # Translate complete raster by x + 100
@@ -132,12 +135,24 @@ julia> GeoArray(rand(5,5,1)) - GeoArray(rand(5,5,1))
 5x5x1 Array{Float64,3} with AffineMap([1.0 0.0; 0.0 1.0], [0.0, 0.0]) and undefined CRS
 ```
 
-### Nodata filling
+One can also warp an array, using GDAL behind the scenes.
+For example, we can vertically transform from the ellipsoid
+to the EGM2008 geoid using EPSG code 3855.
+Note that the underlying PROJ library needs to find the geoidgrids,
+so if they're not available locally, one needs to set `ENV["PROJ_NETWORK"] = "ON"`
+as early as possible, ideally before loading GeoArrays.
+```julia
+ga = GeoArray(zeros((360, 180)))
+bbox!(ga, (min_x=-180, min_y=-90, max_x=180, max_y=90))
+crs!(ga, GeoFormatTypes.EPSG(4979))  # WGS83 in 3D (reference to ellipsoid)
+ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
+```
 
+### Nodata filling
 GeoArrays with missing data can be filled with the [`fill!`](@ref) function.
 
 ```julia
-julia> using GeoEstimation  # or any esimation solver from the GeoStats ecosystem
+julia> using GeoStatsSolvers  # or any estimation solver from the GeoStats ecosystem
 julia> ga = GeoArray(Array{Union{Missing, Float64}}(rand(5, 1)))
 julia> ga.A[2,1] = missing
 [:, :, 1] =
@@ -156,7 +171,6 @@ julia> GeoArrays.fill!(ga, IDW(:band => (neighbors=3,)))  # band is the hardcode
 ```
 
 ### Plotting
-
 Individual bands from a GeoArray can be plotted with the `plot` function. By default the first band is used.
 
 ```julia
@@ -177,7 +191,6 @@ By default the sample size is twice figure size. You can control this factor by
 where higher scalefactor yields higher sizes, up to the original GeoArray size.
 
 ### Subsetting arrays
-
 GeoArrays can be subset by row, column and band using the array subsetting notation, e.g. `ga[100:200, 200:300, 1:2]`.
 
 ```julia
@@ -203,11 +216,9 @@ julia> plot(ga_sub)
 ### Profile
 You can sample the values along a line in a GeoArray with `profile(ga, linestring)`. The linestring can be any geometry that supports [GeoInterface.jl](https://github.com/JuliaGeo/GeoInterface.jl/).
 
-## Reference
-```@autodocs
-Modules = [GeoArrays]
-Order   = [:function, :type]
-```
+
+## Alternatives
+GeoArrays.jl was written to quickly save a geospatial Array to disk. Its functionality mimics `rasterio` in Python. If one requires more features---such as rasterization or zonal stats---which also work on NetCDF files, [Rasters.jl](https://github.com/rafaqz/Rasters.jl/) is a good alternative. Its functionality is more like `(rio)xarray` in Python.
 
 ## Index
 ```@index

src/GeoArrays.jl

@@ -6,6 +6,8 @@ using CoordinateTransformations
 using StaticArrays
 import GeoInterface as GI
 using IterTools: partition
+import DataAPI
+import Extents
 
 include("geoarray.jl")
 include("geoutils.jl")
@@ -15,6 +17,7 @@ include("interpolate.jl")
 include("crs.jl")
 include("operations.jl")
 include("plot.jl")
+include("geointerface.jl")
 
 export GeoArray