rafaqz / Rasters.jl

Raster manipulation for the Julia language
MIT License
197 stars 34 forks source link

`cellsize` should error informatively if there is no CRS #660

Closed rafaqz closed 1 month ago

rafaqz commented 1 month ago

This isn't so nice:

julia> lon = X(-179.875:1:179.875);

julia> lat = Y(-89.875:1:89.875);

julia> ras = Raster(rand(lon, lat));

julia> Rasters.cellsize(ras)
ERROR: MethodError: Cannot `convert` an object of type Nothing to an object of type CoordSys

Closest candidates are:
  convert(::Type{T1}, ::T2) where {T1<:GeoFormat, T2<:T1}
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/UFNTK/src/GeoFormatTypes.jl:93
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84
  convert(::Type{T}, ::AbstractString) where T<:GeoFormat
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/UFNTK/src/GeoFormatTypes.jl:154
  ...

Stacktrace:
 [1] cellsize(dims::Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}})
   @ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/cellsize.jl:54
 [2] cellsize(x::Raster{Float64, 2, Tuple{X{…}, Y{…}}, Tuple{}, Matrix{Float64}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing})
   @ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/cellsize.jl:83
 [3] top-level scope
   @ REPL[70]:1
Some type information was truncated. Use `show(err)` to see complete types.

@tiemvanderdeure

tiemvanderdeure commented 1 month ago

Surely there are a bunch of other functions that should behave similarly? They should probably all throw the same error? Or in some cases warn and assume WGS84?

For example this weird resample just gives a warning

lon = X(-179.875:1:179.875);
lat = Y(-89.875:1:89.875);
ras = Raster(rand(lon, lat));
resample(ras; crs = EPSG(4326))
┌ Warning: You have set a crs to resample to, but the object does not have crs so GDAL will assume it is already in the target crs. Use `newraster = setcrs(raster, somecrs)` to fix this.
└ @ RastersArchGDALExt C:\Users\tsh371\.julia\packages\Rasters\GI5GA\ext\RastersArchGDALExt\resample.jl:60

But reproject(ras; crs = EPSG(4326)) does not even throw a warning and just returns ras.

rafaqz commented 1 month ago

Surely there are more! This is just one that came up for someone. Something consistent would be good.