[ANN] Announcing GeoArrays.jl

evetion · March 8, 2021, 2:50pm

I’d like to announce GeoArrays.jl. To be fair, it’s from 2018, but I think it’s good to have it indexed here.

GeoArrays provides simple geographical raster interaction built on top of ArchGDAL, GDAL and CoordinateTransformations. It follows the design from the excellent rasterio package in Python. That means a GeoArray has an 3D Array (x,y,bands) with data values, an AffineMap to map logical indices to coordinates and a CRS to link the coordinates to the real world.

Reading:

julia> fn = download("https://github.com/yeesian/ArchGDALDatasets/blob/master/data/utmsmall.tif?raw=true")
julia> geoarray = GeoArrays.read(fn)
100x100x1 Array{UInt8,3} with AffineMap([60.0 0.0; 0.0 -60.0], [440720.0, 3.75132e6]) and CRS PROJCS["NAD27 / UTM zone 11N"...

Writing:

julia> ga = GeoArray(rand(100,200))
julia> bbox!(ga, (min_x=2., min_y=51., max_x=5., max_y=54.))  # roughly the Netherlands
julia> epsg!(ga, 4326)  # in WGS84
julia> GeoArrays.write!("test.tif", ga)

An alternative package with a broader scope, such as including netCDFs is GeoData.

juliohm · March 8, 2021, 4:04pm

Similar to the post [ANN] Announcing GeoDataFrames.jl it would be nice to see a GeoArray type implementing the Data trait from Meshes.jl in case there is any interest in doing geostatistics smoothly in that space. I am positive that Meshes.jl needs to be more aware of coordinates tranformations and projections, something that you could certainly contribute with your expertise. We are still brainstorming the possible designs for handling projections of geometries and lat/lon mappings for example.

evetion · March 8, 2021, 6:35pm

Hi Júlio, at the moment GeoArrays implements the old(?) GeoStatsBase for interpolation. To be fair, I haven’t kept up with the latest (breaking) changes, but I love that we can do easy interpolation on the missing values of a GeoArray by using any solver from the GeoStats ecosystem.

In the end, a raster is just a structured mesh or nicely aligned points and it should be easy to conver to such a thing with traits. I’ll look into your Data trait.

juliohm · March 8, 2021, 7:10pm

As per your suggestion, we are migrating a lot of code from GeoStatsBase.jl to other packages such as Meshes.jl. So in theory all you need now is depend on Meshes.jl to implement the trait and then load GeoStats.jl at runtime to do the estimation, simulation etc. Exciting times for these projects. Lots more coming up.

dlb · April 28, 2022, 3:31pm

I’m a newbie and would really appreciate some assistance with reading a very large (250 GB) GeoTIFF file.
Direct load:
GeoArrays.read(“/Users/dB/Downloads/nationalz54ag.tif”)
…kills the Julia session, so try:

ga=GeoArrays.read(“/Users/dB/Downloads/nationalz54ag.tif”, masked=false)
ga.A
GDAL Dataset (Driver: GTiff/GeoTIFF)
*File(s): *

/Users/dB/Downloads/nationalz54ag.tif*

Dataset (width x height): 105002 x 642397 (pixels)
Number of raster bands: 1

[GA_ReadOnly] Band 1 (Gray): 105002 x 642397 (Float32)*

Try to read subset, but it returns nonsense no matter where I look in the raster:

ga.A[1:5,:1:5,:1]
5×5 Matrix{Float32}:

-3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38*
-3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38*
-3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38*
-3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38*
-3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38 -3.40282f38*

Metadata looks sane:

105002x642397x1 ArchGDAL.RasterDataset{Float32, ArchGDAL.IDataset} with AffineMap([5.0 0.0; 0.0 -5.0], [254989.57538994798, 8.9059849999e6]) and CRS PROJCS[“GDA94 / MGA zone 54”,GEOGCS[“GDA94”,DATUM[“Geocentric_Datum_of_Australia_1994”,SPHEROID[“GRS 1980”,6378137,298.257222101,AUTHORITY[“EPSG”,“7019”]],AUTHORITY[“EPSG”,“6283”]],PRIMEM[“Greenwich”,0],UNIT[“degree”,0.0174532925199433,AUTHORITY[“EPSG”,“9122”]]],PROJECTION[“Transverse_Mercator”],PARAMETER[“latitude_of_origin”,0],PARAMETER[“central_meridian”,141],PARAMETER[“scale_factor”,0.9996],PARAMETER[“false_easting”,500000],PARAMETER[“false_northing”,10000000],UNIT[“metre”,1,AUTHORITY[“EPSG”,“9001”]],AXIS[“Easting”,EAST],AXIS[“Northing”,NORTH]]

Any clues?
Thnx dB

juliohm · April 28, 2022, 4:23pm

Take a look at

https://github.com/tlnagy/TiffImages.jl

It is a pure Julia TIFF reader. It should also support GeoTIFF.

After you load your image into an array, you can proceed and work with the other packages. For example, given an image as an array you can georef((img=img,)) to perform geostatistical analysis.

ayushpatnaikgit · April 28, 2022, 4:51pm

Not sure why this is happening, but give ArchGDAL a shot.
https://yeesian.com/ArchGDAL.jl/latest/rasters/

You can do:

using ArchGDAL
img = ArchGDAL.readraster(filepath)

I notice that your image has just 1 band. You can do the following to subset:
Suppose want to crop from:

top_left = [100, 100]
bottom_right = [200, 200]
# These correspond to 2 pairs latitude and longitude on earth.

Then the command will be:

img = img[top_left[2]+1:bottom_right[2], top_left[1]+1:bottom_right[1], 1]

This cropping will not be done in memory, so your Julia session shouldn’t get killed.

I would like to point out that you may need to transpose the image before using it if you end up using ArchGDAL.

I’ve never used a 250 GB image, but I have successfully used this method on images up to 10 GB.

Good luck!

joa-quim · April 28, 2022, 6:18pm

Big GeoTIFF files often come with pyramid layers exactly to avoid having to load the larger array all in memory when that is not needed. But GDAL also allows reading only a sub-region of the file. I’m certain ArchGDAL can do it (just don’t know how) and GMT too, see the gdaltranslate function.

evetion · April 28, 2022, 6:29pm

Try to read subset, but it returns nonsense no matter where I look in the raster:

@dlb Because you’re only loading it lazily, the data is not masked yet. So what you’re seeing are probably the nodata values stored in the raster? Maybe check in the middle of the grid.

You could get the nodatavalue, if it is set for the band with:

const AG = GeoArrays.ArchGDAL
AG.getnodatavalue(AG.getband(ga.A, 1))

It is a pure Julia TIFF reader. It should also support GeoTIFF.

Yes, but without any interpretation of the spatial geo information stored in the TIFF, such as the actual spatial reference and transformation.

Not sure why this is happening, but give ArchGDAL a shot.

GeoArrays actually wraps ArchGDAL, in the example you see that using masked=true returns a ArchGDAL.RasterDataset. So that would give you the same problem.