yeesian / ArchGDAL.jl

A high level API for GDAL - Geospatial Data Abstraction Library
https://yeesian.github.io/ArchGDAL.jl/stable/
Other
137 stars 25 forks source link

GDAL.gdalwarp fails: Point outside of projection domain #373

Closed Rapsodia86 closed 1 year ago

Rapsodia86 commented 1 year ago

Although the issue occurred when I was using GDAL.jl, but since ArchGDAL is built on top of it, I am posting this problem here as well.

https://github.com/JuliaGeo/GDAL.jl/issues/149

UPDATE 2: I reprojected the file using osgeo gdal.Wrap in python 3.8.13, osgeo gdal =3.6.2 And although it showed me the same error: "ERROR 1: Point outside of projection domain", the output was created and looks ok. Any way to have that error suppressed and have the reprojection done?

Using osgeo gdal '3.4.2', python 3.8.13 there is no error but both files have slightly different extents.

UPDATE: Ok, I think the problem is with the source projection ETRS89-extended / LAEA Europe, the newer GDAL version might have a problem. In R terra, I am getting when reading the data:

proj_create_from_database: prime meridian not found
proj_create_from_database: datum not found
proj_create_from_database: ellipsoid not found

Might be related to: https://github.com/r-spatial/stars/issues/588


Hi there, I have just updated my packages in Julia and got an error while running gdalwarp.

options22 = GDAL.gdalwarpappoptionsnew(["-co","COMPRESS=LZW","-srcnodata",string(nodata),"-dstnodata",string(nodata),"-t_srs",pr_m, "-of","GTiff","-tr","100","100",pr_m,  "-r", "nearest"], C_NULL)  
ds_warped22 = GDAL.gdalwarp(fn_out22, Ptr{GDAL.GDALDatasetH}(C_NULL), 1, [vpp2], options22, C_NULL)
ERROR: GDALError (CE_Failure, code 1):
        Point outside of projection domain

Stacktrace:
 [1] gdaljl_errorhandler(class::GDAL.CPLErr, errno::Int32, errmsg::Cstring)
   @ GDAL C:\Users\monikat\.julia\packages\GDAL\wjbw2\src\error.jl:36
 [2] gdalwarp(pszDest::String, hDstDS::Ptr{Ptr{Nothing}}, nSrcCount::Int64, pahSrcDS::Vector{Ptr{Nothing}}, psOptions::Ptr{Nothing}, pbUsageError::Ptr{Nothing})
   @ GDAL C:\Users\monikat\.julia\packages\GDAL\wjbw2\src\libgdal.jl:20498
 [3] top-level scope
   @ Untitled-1:55
 pr_m
"PROJCS[\"Pulkovo 1942(58) / Stereo70\",GEOGCS[\"Pulkovo 1942(58)\",DATUM[\"Pulkovo_1942_58\",SPHEROID[\"Krassowsky 1940\",6378245,298.3,AUTHORITY[\"EPSG\",\"7024\"]],AUTHORITY[\"EPSG\",\"6179\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4179\"]],PROJECTION[\"Oblique_Stereographic\"],PARAMETER[\"latitude_of_origin\",46],PARAMETER[\"central_meridian\",25],PARAMETER[\"scale_factor\",0.99975],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",500000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"3844\"]]"
source projection: GDAL.gdalgetprojectionref(vpp2)
"PROJCS[\"ETRS89-extended / LAEA Europe\",GEOGCS[\"ETRS89\",DATUM[\"European_Terrestrial_Reference_System_1989\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6258\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4258\"]],PROJECTION[\"Lambert_Azimuthal_Equal_Area\"],PARAMETER[\"latitude_of_center\",52],PARAMETER[\"longitude_of_center\",10],PARAMETER[\"false_easting\",4321000],PARAMETER[\"false_northing\",3210000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Northing\",NORTH],AXIS[\"Easting\",EAST],AUTHORITY[\"EPSG\",\"3035\"]]"

Before the update, everything was working fine. now that error. And I do not have the old manifest saved... However, I found under the old Artifact.toml in GDAL_jll was GDAL-vv300.400.100, the new one is v301.600.200 if that helps... Previous GDAL was: 1.3.0

Thanks!

Here are the actual packages installed: Julia v1.7.2 [c9ce4bd3] ArchGDAL v0.9.3 [336ed68f] CSV v0.10.9 [5ae59095] Colors v0.12.10 [717857b8] DSP v0.7.8 [a93c6f00] DataFrames v1.5.0 [31c24e10] Distributions v0.25.86 [add2ef01] GDAL v1.5.1 [c91e804a] Gadfly v1.3.4 [c27321d9] Glob v1.3.0 [cd3eb016] HTTP v1.7.4 [c8e1da08] IterTools v1.4.0 [682c06a0] JSON v0.21.3 [b946abbf] NaNStatistics v0.6.18 [e7bfaba1] NumericalIntegration v0.3.3 [f27b6e38] Polynomials v3.2.5 [438e738f] PyCall v1.95.1 [ce6b1742] RDatasets v0.7.7 [a3a2b9e3] Rasters v0.5.1 [8e980c4a] Shapefile v0.9.1 [2913bbd2] StatsBase v0.33.21 [bd369af6] Tables v1.10.0 [ade2ca70] Dates [10745b16] Statistics

evetion commented 1 year ago

So there are actually two problems:

  1. GDAL (or actually PROJ behind it) makes breaking changes all the time in terms of transformation results. I've been bitten by this before (https://github.com/JuliaGeo/Proj.jl/issues/71), which essentially boils down to, fix your PROJ_jll version or make sure you have tests in place for transformations.
  2. It seems that one can get valid output, even though GDAL throws an error. What happens if you do a try/catch? Do you get something? Otherwise we might need to rethink the error macros in GDAL.jl.
Rapsodia86 commented 1 year ago

Hi, thanks for the answer!

  1. Yeah, breaking changes may really catch you surprised.
  2. We have three situations here: A. Python '3.8.13' Osgeo GDAL '3.6.2': Throws an error: "ERROR 1: Point outside of projection domain" but produces the output file. B. Python '3.8.13' Osgeo GDAL '3.4.2': No error and output produced (small differences compared to output from A) C. Julia '1.7.2.' GDAL.jl '1.5.1': it throws an error and does not produce the output. After throwing the error everything stops. Do you mean basic try/catch? Like this?
    try 
    ds_warped1 = GDAL.gdalwarp(data_out, Ptr{GDAL.GDALDatasetH}(C_NULL), 1, [vpp], options1, C_NULL)   
    catch e
    end

    Well, it catches the error as it should and does not print it (no print function here obviously), and that is all. Unless you do mean some other error-handling approach?

Rapsodia86 commented 1 year ago

"Otherwise we might need to rethink the error macros in GDAL.jl." Anything in that matter? Could that particular error be treated as a warning? Since Osgeo GDAL python does not stop the whole code?

Rapsodia86 commented 1 year ago

Hello, In case of those errors can we have them non-breaking process? Like in case A: A. Python '3.8.13' Osgeo GDAL '3.6.2': Throws an error: "ERROR 1: Point outside of projection domain" but produces the output file.

maxfreu commented 1 year ago

Can you test it again please, now that https://github.com/JuliaGeo/GDAL.jl/pull/153 is merged?

Rapsodia86 commented 1 year ago

Great, I got the result without any error or warning! Thank you!