Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rasters.write for ArchGDAL tiff incorrect blocksize #527

Closed
kongdd opened this issue Sep 28, 2023 · 18 comments · Fixed by #528
Closed

Rasters.write for ArchGDAL tiff incorrect blocksize #527

kongdd opened this issue Sep 28, 2023 · 18 comments · Fixed by #528

Comments

@kongdd
Copy link
Contributor

kongdd commented Sep 28, 2023

using Rasters
using ArchGDAL

using GDAL_jll
function gdal_info(f)
  run(`$(gdalinfo_path()) $f`)
  nothing
end

n = 200
x = X(1:n*7)
y = Y(1:n*4)

A = rand(Int, 7 * n, 4 * n)
ra = Raster(A, dims=(x, y))

case01

fout = "test-terra-03.tif"
@time write(fout, ra; force=true);
# 0.144271 seconds (13.99 k allocations: 34.854 MiB, 7.58% gc time)
gdal_info(fout)
Driver: GTiff/GeoTIFF
Files: test-terra-03.tif
Size is 1400, 800
Origin = (1.000000000000000,801.000000000000000)   
Pixel Size = (1.000000000000000,-1.000000000000000)
Image Structure Metadata:
  COMPRESSION=ZSTD
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       1.000,     801.000)
Lower Left  (   1.0000000,   1.0000000)
Upper Right (    1401.000,     801.000)
Lower Right (    1401.000,       1.000)
Center      (     701.000,     401.000)
Band 1 Block=1400x1 Type=Int64, ColorInterp=Gray   
  NoData Value=-9223372036854775808

The BLOCKYSIZE is always 1. Is this behavior expected?

case02

options = Dict("BLOCKXSIZE"=>"512", "BLOCKYSIZE"=>"512")
@time write(fout, ra; force=true, options); # CO
julia> @time write(fout, ra; force=true, options); # CO
┌ Warning: `missing` cant be written to .tif, missinval for `Int32` of `-2147483648` used instead
└ @ Rasters \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\src\utils.jl:32
(driver, gdaldriver) = ("COG", Driver: COG/Cloud optimized GeoTIFF generator)
options = Dict("COMPRESS" => "ZSTD", "BLOCKXSIZE" => "512", "BLOCKYSIZE" => "512")
ERROR: ArgumentError: Invalid driver creation option(s) detected.
        Please check them carefully with the documentation at https://gdal.org/drivers/raster/index.html.     
        ["BLOCKXSIZE=512", "BLOCKYSIZE=512"]
Stacktrace:
  [1] _gdal_process_options(driver::String, options::Dict{String, String}; _block_template::Raster{Int32, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Matrix{Int32}, Symbol, DimensionalData.Dimensions.LookupArrays.NoMetadata, Int32})
    @ RastersArchGDALExt \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:448
  [2] _gdal_process_options
    @ \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:401 [inlined]       
  [3] _gdal_with_driver(f::RastersArchGDALExt.var"#28#30"{Raster{Int32, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Matrix{Int32}, Symbol, DimensionalData.Dimensions.LookupArrays.NoMetadata, Int32}}, filename::String, driver::String, create_kw::NamedTuple{(:width, :height, :nbands, :dtype), Tuple{Int64, Int64, Int64, DataType}}; options::Dict{String, String}, _block_template::Raster{Int32, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, StepRange{Int64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Matrix{Int32}, Symbol, DimensionalData.Dimensions.LookupArrays.NoMetadata, Int32})
    @ RastersArchGDALExt \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:378
  [4] _gdal_with_driver
    @ \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:373 [inlined]
  [5] #_gdalwrite#27
    @ \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:361 [inlined]
  [6] _gdalwrite
    @ \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:356 [inlined]
  [7] write(filename::String, ::Type{Rasters.GDALsource}, A::Raster{Int32, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, 
UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Matrix{Int32}, Symbol, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing}; force::Bool, verbose::Bool, kw::Base.Pairs{Symbol, Dict{String, String}, Tuple{Symbol}, NamedTuple{(:options,), Tuple{Dict{String, String}}}})
    @ RastersArchGDALExt \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:44
  [8] write
    @ \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\ext\RastersArchGDALExt\gdal_source.jl:35 [inlined]
  [9] write(filename::String, A::Raster{Int32, 2, Tuple{X{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}, Y{DimensionalData.Dimensions.LookupArrays.Sampled{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Int64}, DimensionalData.Dimensions.LookupArrays.Points, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Tuple{}, Matrix{Int32}, Symbol, DimensionalData.Dimensions.LookupArrays.NoMetadata, Missing}; source::Type, kw::Base.Pairs{Symbol, Any, Tuple{Symbol, Symbol}, NamedTuple{(:force, :options), Tuple{Bool, Dict{String, String}}}})
    @ Rasters \\kong-nas\CMIP6\GitHub\jl-spatial\Rasters.jl\src\write.jl:11
 [10] top-level scope
    @ .\timing.jl:273 [inlined]
 [11] top-level scope
    @ .\REPL[33]:0

package version

  • Rasters v0.8.4
  • ArchGDAL v0.10.1
@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

The tiff file doc: https://gdal.org/drivers/raster/gtiff.html

image

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Re-wrapping Rasters.jl and calling it Terra.jl is pretty weird... why do you do things like this all the time?

Also, please post a full stack trace so I an actually find where that error is throw instead of scrolling around looking for it.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

Revised. This issue has no relationship with Terra.jl.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Haha but why does Terra.jl exist at all is the question.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

You still havent included a full stack trace either.

If you insist in continuing to wrap my work as random new packages, you can at least make your issues easy to read.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

The above script has been updated, #527 (comment).

This issue is not related to Terra.jl. click the above link, you will see the revised version.
I only used the gdal_info function to print tiff info. Now, I copied it to here.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

Haha but why does Terra.jl exist at all is the question.

In the R user, I want to use Rasters.jl in my way, that you and others may not accept.
Push my function to Rasters.jl, you might refuse and this workflow is low efficient. Terra.jl is a solution for my case. If you don't like, I will not publish it to Julia officical site.
I host it on the github, becasue only public repo have sufficient github action hours.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

@kongdd you really need to explain why you keep making these new packages that add very little value.

It looks like a scam to make it look like our work was done by you.

You code is fine and you could be a respected part of an open source community, but you currently have a reputation as a scammer trying to get credit for other peoples work.

Edit: we were writing at the same time, you did explain some of that.

But if the package is just for you, why call it Terra.jl

Calling it Terra.jl is looking for publicity and clicks - that suggests it is not just a package for you to use.

I rarely refuse PRs, and functions like your gdal_info would actually be useful in Rasters.jl

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

I earned github stars from my R, GEE packages.
I need Julia to do some computation for my research. I revised some function to make it more convenient for me. I didn't mean to get credit from other peoples work.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Well, you need to stop using names like Stars.jl, Terra.jl, and hosting other peoples packages in your organisation.

To someone who doesn't understand github, it may look like all this is your work. And that may be good for your academic career.

That is what many people in the julia community suspect you are trying to do.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

I will delete Stars.jl after all my peasonal function has been migrated.
But for Terra.jl, I really don't understand the community of Julia. Most scripts are written by myself. Indeed, I import the function Rasters.jl. I can declare it in the readme file.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Absolutely declare that it's a wrapper of Rasters.jl in the Readme.

This is an example of a good use of Rasters.jl:
https://github.com/JoshuaBillson/RemoteSensingToolbox.jl

Its mentioned immediately in the readme, it adds a heap of field specific functions for people to use in remote sensing, and has a name that lets people know that that is what it does. And the author contributes code to DimensionalData.jl

Calling your package Terra.jl is different. Its saying to google search that its like terra in R. Which is not true.

Its also saying it is a direct alternative to Rasters.jl, so if people use it you are essentially forking the ecosystem to have two packages with similar components and lots of duplicated work.

JuliaGeo is actually just a few people trying to put all these tools together, and we actually need help and more packages that integrate with what we have built, rather than duplicating and confusing it.

I have really appreciated your pull request in the past and would love to get them in future, or to see new packages that add cool field specific tools to Rasters.jl like RemoteSensingToolbox does.

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Finally, getting back to the issue...

The problem is GDAL recently switched the default backend for .tif files to "COG", which uses "BLOCKSIZE" rather than having separate X/Y values.

This causes a bunch of problems because "COG" is not a good default for lazy writing, so we have to use "GeoTIFF" anyway.

So it seems the problem is your options are correct but we are still checking them against the allowed arguments for "COG" instead of "GeoTIFF". It should be easy to fix.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

Declaration is added in the README file. Further suggestion is welcome.
Stars.jl is removed.

Everyone want his tools has a nice name, I currently have no idea how to rename Terra.jl.

image

@rafaqz
Copy link
Owner

rafaqz commented Sep 28, 2023

Thanks, but still... using names of common packages will make people either confused or suspicious,
unless you are wrapping a C binary like GDAL.jl does. Then the name makes it clear what the package is.

Rasters.jl has a name like Rs "raster" because that's the common name in english of the object we all work with - a matrix of pixels is a "raster". rasterio in python is similar.

But terra is a nice name that means "earth" that someone thought of for that specific R package. Its not a common word for the object or functions like raster is.

So using terra as a name for your package really sounds like your are copying the package or associated with them somehow, or that its a julia wrapper for the R code. It will also get google search results for people wanting that package specifically.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

After revised, this issue has no relationship with 'Terra.jl'. Although I can feel your rudeness, I admit that you are a great man. You contribute a lot to the community.

While, I am lazy. I just want to solve my own problems quickly. Why I should push PRs to Rasters.jl? (1) For me, it is inefficient. After 2-3 days, or 2-3 weeks, after this PR accept, then turn back to my problem? (2) Isn't that a form of intellectual plagiarism? I can't found where is my previous committed script. (3) In my repo, I take control of everything in my favorite. For example, I want to use Raster in the way of rast(A, bbox). But I am not sure whether you will accept it or not. Maybe 2-3 days passed, and then refused.

About the misleading, I admit that I should acknowledge your package, your work and states what this package have. This is my bad. Sorry for my laziness and non-fluent English. But I have the freedom how to name my package. I can change it if there are good candidates.
But it is a personal package. I can give clear declaration. I don't mind about the stars and github reputation. If you offer some sponsorship for PRO GITHUB, I am glad to turn my package into private.

Most R users only know Julia is japanese famous star. Who care about it?

ps. If you can speak in Chinese, I can reply more quickly with more details.

@kongdd
Copy link
Contributor Author

kongdd commented Sep 28, 2023

Anyway, thanks for your effort and time for solving this issue. \heart

One line conclusion: jl-spatial renamed to jl-pkgs, Stars.jl delete, Terra.jl give clear declaration of its dependency of Rasters.

@rafaqz
Copy link
Owner

rafaqz commented Sep 29, 2023

But I am not sure whether you will accept it or not. Maybe 2-3 days passed, and then refused

I review pretty fast here and try to be very open and helpful wherever possible. But, what you are describing is just what I have to do constanly anyway:

2023-09-29-110859_368x227

I certainly don't have 93 packages or projects on the go. Half of those are other peoples packages that I have to contribute to to make things work here and in my other packages.

Often no one reviews for a week or more, because they are busy. That's just what its like working collectively.

Isn't that a form of intellectual plagiarism?

Commits always have your name on them - you are forever recognized as a contributor here:

2023-09-29-111200_597x356

You can absolutely work alone however you like. I also have my own projects I don't register. The problem is if you try to publicize the work, like when you tried to register Stars.jl when it was not your work. That makes us worried you will do it again.

Imagine if I took one of your papers, took your name off it, and published it in a different journal. Would that be ok with you?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants