evetion / GeoArrays.jl

Simple geographical raster interaction built on top of ArchGDAL, GDAL and CoordinateTransformations
https://www.evetion.nl/GeoArrays.jl/dev/
MIT License
51 stars 13 forks source link

Local test failure #161

Closed hustf closed 4 months ago

hustf commented 4 months ago

It seems a geoarray is all zeros when it shouldn't. Why does it happen on my computer and not in the automated tests? My context is shown at the end.

Reported as when running #160 locally. I haven't pulled the last two commits, just the initial [3cfcd27].

]test
...
warp: Test Failed at C:\Users\f\.julia\dev\GeoArrays\test\test_operations.jl:60
  Expression: sum(ga2) != 0
   Evaluated: 0.0 != 0
...
  operations                           |   23     1            24   3.4s
    dimension error test               |    4                   4   0.9s
    affine error tests                 |    4                   4   0.0s
    crs error test                     |    4                   4   0.0s
    operations                         |    9                   9   1.0s
    warp                               |    2     1             3   1.4s
  getdriver                            |                     None   0.0s
ERROR: LoadError: Some tests did not pass: 356 passed, 1 failed, 1 errored, 0 broken.
in expression starting at C:\Users\f\.julia\dev\GeoArrays\test\runtests.jl:8
ERROR: Package GeoArrays errored during testing

I then tried running this locally. using GeoArrays, etc.:

julia> pwd()
"c:\\Users\\f\\.julia\\dev\\GeoArrays\\test"

julia> begin
        ga = GeoArray(zeros((360, 180)))
        bbox!(ga, (min_x=-180, min_y=-90, max_x=180, max_y=90))
        crs!(ga, GeoFormatTypes.EPSG(9754))
        ga2 = GeoArrays.warp(ga, Dict("t_srs" => "EPSG:4326+3855"))
end
360x180x1 Array{Float64, 3} with AffineMap([1.0 0.0; 0.0 -1.0], [-180.0, 90.0]) and CRS COMPD_CS["WGS 84 + EGM2008 height",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]],VERT_CS["EGM2008 height",VERT_DATUM["EGM2008 geoid",2005,AUTHORITY["EPSG","1027"]],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Gravity-related height",UP],AUTHORITY["EPSG","3855"]]]

julia> sum(ga2)  # The failed test: This is the wrong answer!
0.0

julia> sum(abs.(ga2.A)) # ...and not by coincidence.
0.0

julia> VERSION
v"1.10.2"

(GeoArrays) pkg> status
Project GeoArrays v0.8.5
Status `C:\Users\f\.julia\dev\GeoArrays\Project.toml`
⌃ [c9ce4bd3] ArchGDAL v0.10.3 ⚲
  [150eb455] CoordinateTransformations v0.6.3
  [9a962f9c] DataAPI v1.16.0
  [411431e0] Extents v0.1.2
  [68eda718] GeoFormatTypes v0.4.2
  [cf35fbd7] GeoInterface v1.3.4
  [c8e1da08] IterTools v1.10.0
  [aea7be01] PrecompileTools v1.2.1
  [3cdcf5f2] RecipesBase v1.3.4
  [90137ffa] StaticArrays v1.9.4
⌃ [a7073274] GDAL_jll v301.900.0+0
⌃ [58948b4f] PROJ_jll v901.300.0+1

julia> Sys.KERNEL  # This on Windows 64-bit
:NT

julia> Sys.CPU_THREADS  # Julia's running with 12 threads. Running with 1 thread: same answer.
12

julia> function get_system_info(command)
           return read(`cmd /c $command`, String)
       end;

julia> get_system_info("wmic cpu get name, NumberOfCores, NumberOfLogicalProcessors /format:list") |> println

Name=AMD Ryzen 5 3600X 6-Core Processor
NumberOfCores=6
NumberOfLogicalProcessors=12
evetion commented 4 months ago

Thanks for making an issue! I'm sorry you hit this issue, it's caused because Proj silently fails to find the vertical transformation and thus does nothing. There's an upstream issue to change this behaviour.

Our CI sets an ENV variable, so Proj downloads the missing data on the fly https://github.com/evetion/GeoArrays.jl/blob/master/.github%2Fworkflows%2FCI.yml#L46

hustf commented 4 months ago

Ok, upstream issues are inevitable! I hope you get as nice a response as I got!

By the way, I did set ENV["PROJ_NETWORK"] = "ON" as in the runtests.jl, but that's non-related I suppose.

If you haven't merged yet, I'll review the PR though I'm not qualified regarding the rest of the package....

evetion commented 4 months ago

By the way, I did set ENV["PROJ_NETWORK"] = "ON" as in the runtests.jl, but that's non-related I suppose.

I fear it needs to be set before you load (either) GeoArrays/ArchGDAL/GDAL/Proj, as it is checked at the initialization of the library, which happens as soon as you do using. You can use this shortcut to enable network access: GeoArrays.ArchGDAL.GDAL.osrsetprojenablenetwork(true). I'm looking for a way to enable it as soon as one requests a vertical transformation.