Merge branch 'as/fixrasterize' into as/makie_and_cf

Project.toml

@@ -82,6 +82,7 @@ ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CFTime = "179af706-886a-5703-950a-314cd64e0468"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
@@ -94,4 +95,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "RasterDataSources", "SafeTestsets", "Shapefile", "Statistics", "Test", "ZarrDatasets"]

docs/src/index.md

@@ -36,17 +36,11 @@ features:
     Rasters provides a standardised interface that allows many source data types to
     be used with identical syntax.
 
-    - Scripts and packages building on Rasters.jl can treat `Raster`,
-    `RasterStack`, and `RasterSeries` as black boxes.
-    - The data could hold GeoTiff or NetCDF files, `Array`s in memory or
-        `CuArray`s on the GPU - they will all behave in the same way.
-    - `RasterStack` can be backed by a Netcdf or HDF5 file, or a `NamedTuple` of
-        `Raster` holding `.tif` files, or all `Raster` in memory.
+    - Scripts and packages building on Rasters.jl can treat `Raster`, `RasterStack`, and `RasterSeries` as black boxes.
+    - The data could hold GeoTiff or NetCDF files, `Array`s in memory or `CuArray`s on the GPU - they will all behave in the same way.
+    - `RasterStack` can be backed by a Netcdf or HDF5 file, or a `NamedTuple` of `Raster` holding `.tif` files, or all `Raster` in memory.
     - Users do not have to deal with the specifics of spatial file types.
-    - `Projected` lookups with Cylindrical projections can by indexed using other Cylindrical projections
-    by setting the `mappedcrs` keyword on construction. You don't need to know the underlying
-    projection, the conversion is handled automatically. This means lat/lon
-    `EPSG(4326)` can be used seamlessly if you need that.
+    - `Projected` lookups with Cylindrical projections can by indexed using other Cylindrical projections by setting the `mappedcrs` keyword on construction. You don't need to know the underlying projection, the conversion is handled automatically. This means lat/lon `EPSG(4326)` can be used seamlessly if you need that.
 
 ## Installation
 
@@ -104,4 +98,4 @@ To make an issue we can fix quickly (or at all) there are three key steps:
 
 Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.
 
-:::
+:::

ext/RastersArchGDALExt/resample.jl

@@ -55,7 +55,16 @@ function resample(A::RasterStackOrArray;
             # get crs from `to` or `x` if none was passed in
             isnothing(Rasters.crs(to)) ? Rasters.crs(A) : Rasters.crs(to)
         end
-    else
+    else # issomething(crs)
+        if crs isa String
+            error("""
+                Strings as CRS aren't yet supported.  
+                Please pass a `GeoFormatTypes.jl` CRS format, like `ESRIWellKnownText`, `ProjString`, or similar.
+
+                You can find out more about GeoFormatTypes.jl at https://juliageo.org/GeoFormatTypes.jl/stable/.
+                """
+            )
+        end
         crs
     end
     if !isnothing(crs)

src/Rasters.jl

@@ -145,17 +145,4 @@ include("sources/grd.jl")
 include("sources/commondatamodel.jl")
 include("extensions.jl")
 
-# Compatibility with pre-1.9 julia
-function __init__()
-    @static if !isdefined(Base, :get_extension)
-        @require ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" include("../ext/RastersArchGDALExt/RastersArchGDALExt.jl")
-        @require CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" include("../ext/RastersCoordinateTransformationsExt/RastersCoordinateTransformationsExt.jl")
-        @require HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" include("../ext/RastersHDF5Ext/RastersHDF5Ext.jl")
-        @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("../ext/RastersMakieExt/RastersMakieExt.jl")
-        @require NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab" include("../ext/RastersNCDatasetsExt/RastersNCDatasetsExt.jl")
-        @require GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc" include("../ext/RastersGRIBDatasetsExt/RastersGRIBDatasetsExt.jl")
-        @require RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36" include("../ext/RastersRasterDataSourcesExt/RastersRasterDataSourcesExt.jl")
-    end
-end
-
 end

src/methods/zonal.jl

@@ -19,7 +19,10 @@ These can be used when `of` is or contains (a) GeoInterface.jl compatible object
     where the line `:touches` the pixel, or that are completely `:inside` inside the polygon.
     The default is `:center`.
 - `progress`: show a progress bar, `true` by default, `false` to hide..
-
+- `skipmissing`: wether to apply `f` to the result of `skipmissing(A)` or not. If `true`
+    `f` will be passed an iterator over the values, which loses all spatial information.
+    if `false` `f` will be passes a masked `Raster` or `RasterStack`, and will be responsible
+    for handling missing values itself. The default value is `true`.
 
 # Example
 
@@ -71,15 +74,18 @@ insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Z
 """
 zonal(f, x::RasterStackOrArray; of, kw...) = _zonal(f, x, of; kw...)
 
-_zonal(f, x::RasterStackOrArray, of::RasterStackOrArray) = _zonal(f, x, Extents.extent(of))
-_zonal(f, x::RasterStackOrArray, of::DimTuple) = _zonal(f, x, Extents.extent(of))
+_zonal(f, x::RasterStackOrArray, of::RasterStackOrArray; kw...) = 
+    _zonal(f, x, Extents.extent(of); kw...)
+_zonal(f, x::RasterStackOrArray, of::DimTuple; kw...) = 
+    _zonal(f, x, Extents.extent(of); kw...)
 # We don't need to `mask` with an extent, it's square so `crop` will do enough.
-_zonal(f, x::Raster, of::Extents.Extent) = f(skipmissing(crop(x; to=of, touches=true)))
-function _zonal(f, x::RasterStack, ext::Extents.Extent)
+_zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
+    _maybe_skipmissing_call(f, crop(x; to=of, touches=true), skipmissing)
+function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true)
     cropped = crop(x; to=ext, touches=true)
     prod(size(cropped)) > 0 || return missing
     return map(cropped) do A
-        f(skipmissing(A))
+        _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
 # Otherwise of is a geom, table or vector
@@ -89,22 +95,28 @@ _zonal(f, x, ::GI.AbstractFeatureCollectionTrait, fc; kw...) =
     _zonal(f, x, nothing, fc; kw...)
 _zonal(f, x::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...) =
     _zonal(f, x, GI.geometry(feature); kw...)
-function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; kw...)
+function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; 
+    skipmissing=true, kw...
+)
     cropped = crop(x; to=geom, touches=true)
     prod(size(cropped)) > 0 || return missing
     masked = mask(cropped; with=geom, kw...)
-    return f(skipmissing(masked))
+    return _maybe_skipmissing_call(f, masked, skipmissing)
 end
-function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom; kw...)
+function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom; 
+    skipmissing=true, kw...
+)
     cropped = crop(st; to=geom, touches=true)
     prod(size(cropped)) > 0 || return map(_ -> missing, st)
     masked = mask(cropped; with=geom, kw...)
     return map(masked) do A
         prod(size(A)) > 0 || return missing
-        f(skipmissing(A))
+        _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
-function _zonal(f, x::RasterStackOrArray, ::Nothing, data; progress=true, threaded=true, geometrycolumn=nothing, kw...)
+function _zonal(f, x::RasterStackOrArray, ::Nothing, data; 
+    progress=true, threaded=true, geometrycolumn=nothing, kw...
+)
     geoms = _get_geometries(data, geometrycolumn)
     n = length(geoms)
     n == 0 && return []
@@ -136,3 +148,5 @@ function _alloc_zonal(f, x, geoms, n; kw...)
     zs[n_missing + 1] = z1
     return zs, n_missing + 1
 end
+
+_maybe_skipmissing_call(f, A, sm) = sm ? f(skipmissing(A)) : f(A)

src/series.jl

@@ -24,6 +24,7 @@ abstract type AbstractRasterSeries{T,N,D,A} <: AbstractDimArray{T,N,D,A} end
 DD.metadata(A::AbstractRasterSeries) = NoMetadata()
 DD.name(A::AbstractRasterSeries) = NoName()
 DD.label(A::AbstractRasterSeries) = ""
+isdisk(A::AbstractRasterSeries) = any(isdisk, A)
 
 """
     modify(f, series::AbstractRasterSeries)

src/stack.jl

@@ -34,7 +34,7 @@ missingval(s::AbstractRasterStack, name::Symbol) = _singlemissingval(missingval(
 filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = map(s -> filename(s), stack)
 filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = filename(parent(stack))
 
-isdisk(st::AbstractRasterStack) = isdisk(layers(st, 1))
+isdisk(st::AbstractRasterStack) = any(isdisk, layers(st))
 
 setcrs(x::AbstractRasterStack, crs) = set(x, setcrs(dims(x), crs)...)
 setmappedcrs(x::AbstractRasterStack, mappedcrs) = set(x, setmappedcrs(dims(x), mappedcrs)...)

test/methods.jl

@@ -300,6 +300,17 @@ end
         zonal(sum, st; of=dims(st)) == 
         zonal(sum, st; of=Extents.extent(st)) == 
         sum(st)
+
+    @testset "skipmissing" begin
+        a = Array{Union{Missing,Int}}(undef, 26, 31)
+        a .= (1:26) * (1:31)'
+        a[1:10, 3:10] .= missing
+        rast = Raster(a, (X(-20:5), Y(0:30)))
+        @test zonal(sum, rast; of=polygon, skipmissing=false) === missing
+        @test zonal(sum, rast; of=polygon, skipmissing=true) isa Int
+        @test !zonal(x -> x isa Raster, rast; of=polygon, skipmissing=true)
+        @test zonal(x -> x isa Raster, rast; of=polygon, skipmissing=false)
+    end
 end
 
 @testset "zonal return missing" begin

test/rasterize.jl

@@ -1,4 +1,5 @@
 using Rasters, Test, ArchGDAL, ArchGDAL.GDAL, Dates, Statistics, DataFrames, Extents, Shapefile, GeometryBasics
+import GeoDataFrames
 import GeoInterface as GI
 using Rasters.Lookups, Rasters.Dimensions 
 using Rasters: bounds
@@ -510,3 +511,18 @@ end
     # Too slow and unreliable to test in CI, but it warns and uses one thread given 32gb of RAM: 
     # coverage(union, shphandle.shapes; threaded=true, res=1, scale=1000)
 end
+
+@testset "`geometrycolumn` kwarg and detection works" begin
+    # Replicate pointtable
+    fancy_table = deepcopy(pointdf)
+    fancy_table.geom = pointdf.geometry
+    select!(fancy_table, Not(:geometry))
+    # Test that rasterization works with provided geometry column
+    # Just test that it works and does not warn.
+    @test_nowarn rasterize(last, fancy_table; to = A1, geometrycolumn = :geom, fill = 1)
+    # Now add GeoDataFrames blessed metadata keys
+    DataFrames.metadata!(fancy_table, "GEOINTERFACE:geometrycolumns", (:geom,); style = :note)
+    # Test that we don't have to provide the geometry column explicitly
+    @test_nowarn rasterize(last, fancy_table; to = A1, fill = 1)
+    @test replace_missing(rasterize(last, pointtable; to = A1, fill = 1), 0) == replace_missing(rasterize(last, fancy_table; to = A1, fill = 1), 0) # sanity check
+end

test/series.jl

@@ -122,5 +122,8 @@ end
         @test all(Rasters.filename.(series) .== filenames)
         first_dims = dims(first(series))
         @test all(dims(r) == first_dims for r in series)
-    end
+        @test Rasters.isdisk(series)
+        @test !Rasters.isdisk(read(series))
+        @test Rasters.isdisk(Rasters.combine(series))
+        end
 end

test/sources/gdal.jl

@@ -615,6 +615,7 @@ end
 
     @testset "lazy" begin
         gdalstack_lazy = RasterStack((a=gdalpath, b=gdalpath); lazy=true)
+        @test Rasters.isdisk(gdalstack_lazy)
         @test Rasters.isdisk(gdalstack_lazy.a)
         @test Rasters.isdisk(gdalstack_lazy.b)
         gdalstack_read = read(gdalstack_lazy)