# Using ArchGDAL

**URL:** <https://discourse.julialang.org/t/using-archgdal/46006>\
**Category:** New to Julia\
**Tags:** question, geo\
**Created:** [September 3, 2020, 1:22pm UTC](https://discourse.julialang.org/t/using-archgdal/46006 "2020-09-03T13:22:27Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [September 3, 2020, 1:22pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/1 "2020-09-03T13:22:27Z")

</div>

Hi, I wanted,

to give Julia a chance because it said to be rather fast. And some parts of my Python code needs hours (like 3 to 5) and that just to much, to be productively useful. Especial when you still have to debug the code.

I want to do some calculation on some masked Sentinel2 tiles and use ArchGDAL.jl for that.

The main code: is opening two raster files and (one with real data, the other as indicies of the reagions. And than I want to compute the mode of each region and than save the data of the mode into an new file.

The code that is currently running (but I wouldn’t mind, when it be faster) is:

```julia
using StatsBase
using ArchGDAL; const AG = ArchGDAL

class_raster = Array{Int64}(ArchGDAL.read(AG.read(prediction_path), 1))
class_raster[class_raster .== 55537] .= 0
class_raster[class_raster .== -9999] .= 0
Array{UInt8}(class_raster)
indices = Array{Int64}(ArchGDAL.read(AG.read(atkis_area_path), 1))

unique_indices = unique(indices)
for _index in unique_indices
       
    coordinates = findall(indices .== _index)
    mode_val = mode(class_raster[findall(index_raster .== _index)])
    for x in coordinates
        local_classes[x] = mode_val
    end
end

```

But when I want to save the file for instance with:

```julia
dataset = AG.read(prediction_path)
ref = AG.getproj(dataset)
geotransform = AG.getgeotransform(dataset)

raster = AG.unsafe_create(
        "test.tif",
        AG.getdriver("GTiff"),
        width = ArchGDAL.width(dataset),
        height = ArchGDAL.height(dataset),
        nbands = 1,
        dtype = UInt8
    )
AG.setgeotransform!(raster, geotransform)
AG.setproj!(raster, ref)

AG.write!(
    raster,
    local_classes, # image to "burn" into the raster
    1, # update band 1
    )

```

I still get an error when calling unsave\_create:

> > MethodError: no method matching unsafe\_create(::String, ::ArchGDAL.Driver; width=7328, height=5860, nbands=1, dtype=UInt8)  
> > Closest candidates are:  
> > unsafe\_create(::AbstractString; driver, width, height, nbands, dtype, options) at C:\Users\TITAN.julia\packages\ArchGDAL\s5Bwr\src\dataset.jl:191
> 
> Stacktrace:  
> [1] top-level scope at In[28]:1  
> [2] include\_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

I don’t know why! Didn’t I add every parameter needed? Is it a problem with the types? Additionally I totally puzzled when to use AG and when ArchGDAL.

I hope you could help me out. Thanks alot!

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 3, 2020, 1:34pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/2 "2020-09-03T13:34:08Z")

</div>

If you’re loading non-rotated regular tif (or other kinds of raster) files, [GeoData.jl](https://github.com/rafaqz/GeoData.jl) wraps ArchGDAL.jl, and can do this in a lot less code:

```julia
using GeoData, ArchGDAL
raster = GDALarray(prediction_path) |> GeoArray
indices = GDALarray(atkis_area_path) |> GeoArray

# Do your things

# Then save
write("output.tif", GDALarray, results)

```

Its a newly published packages so there might be a few bugs, but I use it every day so it works OK.

All the tif metadata is maintained through any math, array slicing or broadcasts. You can just `write` the results array as a tif and it should work. Make an issue on GeoData.jl if it doesn’t.

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [September 3, 2020, 3:10pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/3 "2020-09-03T15:10:13Z")

</div>

> [@Ordens\_Ritter](#):
>
> Additionally I totally puzzled when to use AG and when ArchGDAL.

You can do either as `const AG = ArchGDAL` defines those to be the same. But really you should then just use `AG`, otherwise what’s the point of that definition.

---

<div class="post-metadata">

**Author:** ![visr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/visr/32/17204_2.png) [@visr](https://discourse.julialang.org/u/visr)\
**Post date:** [September 3, 2020, 4:27pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/4 "2020-09-03T16:27:01Z")

</div>

> [@Ordens\_Ritter](#):
>
> MethodError: no method matching unsafe\_create(::String, ::ArchGDAL.Driver; width=7328, height=5860, nbands=1, dtype=UInt8)  
> Closest candidates are:  
> unsafe\_create(::AbstractString; driver, width, height, nbands, dtype, options) at C:\Users\TITAN.julia\packages\ArchGDAL\s5Bwr\src\dataset.jl:191

If you read this very carefully it is saying that you are passing the driver as a positional argument, whereas it is expecting it as a keyword argument, behind the `;`. Where did you base this code on? Would be good to know if something in the docs needs to be updated. There is also no need to use `unsafe_create` anymore, you can directly use `create`, also without using a do-block, see the two examples here:

[https://yeesian.com/ArchGDAL.jl/stable/reference/#ArchGDAL.create-Tuple{AbstractString}](https://yeesian.com/ArchGDAL.jl/stable/reference/#ArchGDAL.create-Tuple%7BAbstractString%7D)

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [September 8, 2020, 4:54pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/5 "2020-09-08T16:54:03Z")

</div>

Hi,  
thanks at @Raf. Your answer has been very fast. Still I could not really apply this method, because I would have to write with an GeoArray… and It didn’t know how to create those… so I installed GeoArrays… and it seamed like using the packages GeoData, ArchGDAL and GeoArray lead to many problems.

For the usage of ArchGDAl’s copy: I neither saw the difference, of the ; and the , now would I have known, what it meaned aktough I read I Julia book so far.

My code currently looks like:

> using ArchGDAL; const AG = ArchGDAL
> 
> prediction\_path = “E:\SeH\8\_TrainingData\output\Python\results\RDF\_Level3\_full\_Prediction\_DLM50\_20m\_manual\_v2\_41features.tif”  
> atkis\_area\_path = “E:\SeH\COP4EE-2\_BB\_indexes.tif”
> 
> raster = Array{Int64}(ArchGDAL.read(AG.read(prediction\_path), 1))  
> raster[raster .== 55537] .= 0  
> raster[raster .== -9999] .= 0  
> raster = Array{UInt8}(raster)  
> indices = Array{Int64}(ArchGDAL.read(AG.read(atkis\_area\_path), 1))
> 
> unique\_indices = unique(indices)
> 
> Threads.@threads for \_index in unique\_indices
> 
> ```
> coordinates = findall(indices .== _index)
> raster[coordinates] .= StatsBase.mode(raster[coordinates])
> 
> ```
> 
> end
> 
> dataset = AG.read(prediction\_path)  
> ref = AG.getproj(dataset)  
> geotransform = AG.getgeotransform(dataset)
> 
> raster\_ds = AG.create(  
> “test.tif”,  
> driver = AG.getdriver(“GTiff”),  
> width = ArchGDAL.width(dataset),  
> height = ArchGDAL.height(dataset),  
> nbands = 1,  
> dtype = UInt8  
> )  
> AG.setgeotransform!(raster\_ds, geotransform)  
> AG.setproj!(raster\_ds, ref)
> 
> AG.write!(  
> raster\_ds,  
> raster, # image to “burn” into the raster  
> 1, # update band 1  
> )

and it’s really faster than python, because I am usage serval threads for the loop.

The only problem I currently have is, that the file that is created seams to be broken. ☹

Does anyone has a clue why?

Thanks alot!

P.S.:

The Output is:

> GDAL Dataset (Driver: GTiff/GeoTIFF)  
> File(s):  
> test.tif
> 
> Dataset (width x height): 7328 x 5860 (pixels)  
> Number of raster bands: 1  
> [GA\_Update] Band 1 (Gray): 7328 x 5860 (UInt8)

But the file is empty.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 9, 2020, 2:06am UTC](https://discourse.julialang.org/t/using-archgdal/46006/6 "2020-09-09T02:06:23Z")

</div>

That script I posted makes a `GeoData.GeoArray`

Thats what this does:

```julia
using GeoData, ArchGDAL
raster = GDALarray(prediction_path) |> GeoArray

```

When prediction\_path is some `String` path to a file GDAL can open.

You can save it again with:

```julia
write("output.tif", GDALarray, raster)

```

Or as a Netcdf:

```julia
using NCDatasets
write("output.ncd", NCDarray, raster[Band(1)])

```

The point is there isn’t anything you have to do.

You can also operate on it in between loading and saving:

```julia
raster = GDALarray(prediction_path; usercrs=EPSG(4326)) |> GeoArray
# Slice a subsection and multiply it by 100
raster = raster[Lat(Between(20, 20), Lon(Between(90, 180))] .* 100
# Save the new raster
write("output.tif", GDALarray, raster)

```

That just works you don’t have to do anything. Most julia array math and methods will just work  
as if it’s a normal array, but it will always remain a GeoArray you can plot and save to disk.

It also doesn’t matter what the underlying CRS is you can just index with `Between` using EPSG:4326 lat/lon coords and it will work if you have set the `usercrs` as above.

You can also plot it:

```julia
using Plots
raster |> plot

```

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [September 14, 2020, 8:48am UTC](https://discourse.julialang.org/t/using-archgdal/46006/8 "2020-09-14T08:48:41Z")

</div>

@Raf Thanks, for the help. The original code idea is working now. Although, I hoped it would be faster even when using Multithreading. Is there a way to change the datatype of the attribute data in the GeoArray?

The new problem is how to add or remove bands to an GeoArray, or can I merge two GeoArray, so that these combine to an GeoArray of two bands?

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 14, 2020, 9:29am UTC](https://discourse.julialang.org/t/using-archgdal/46006/9 "2020-09-14T09:29:53Z")

</div>

What is not fast, loading the file or indexing in to it? Maybe post some code here that is slow and we can work through it (even a link to a file too).

The easiest way to change the array element type is to just broadcast:

```julia
A = Float32.(A)

```

Which will give you a new `GeoArray` holding `Float32`. If that is what you mean?

Joining two arrays in the Band (or any) dimension you can do:

```julia
C = cat(A, B; dims=Band)

```

But they might both be band 1, so the joined dim index will be [1, 1], which you probably don’t want. I’m going to fix `cat` a little to make this easier, but for now you can do:

```julia
C = cat(A, B; dims=Band)
C = setdims(C, Band(1:2, NoIndex()))

```

Although there is a better way of doing this coming with [set](https://github.com/rafaqz/DimensionalData.jl/pull/157)

Edit: It seems to save fine even if the bands are `[1, 1]`, and will load again as `1:2`, so it might not matter.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [September 14, 2020, 10:52am UTC](https://discourse.julialang.org/t/using-archgdal/46006/10 "2020-09-14T10:52:36Z")

</div>

@Raf Thanks for the Tip with the broadcast.

The current speed is limited by the indexing. Especially because it’s over different indices.

It takes about 30 min for this sentinel 2 tile, jupyter using 4 threads on a i7-6700K at 4 GHz.

I splitted the problem in 2 functions:

```julia
function apply_multiband_index_filter(prediction_path::String, index_path::String)
    raster = GDALarray(prediction_path) |> GeoArray
    raster[raster .== 55537] .= 0
    raster[raster .== -9999] .= 0
    raster = UInt8.(raster)
        
    indices = GDALarray(index_path) |> GeoArray
    indices[indices .== 55537] .= 0
    indices[indices .== -9999] .= 0
    indices = Array{Int64}(indices)
    
    county_indices = indices[:, :, 1]
    atkis_indices = indices[:, :, 2]
    invekos_indices = indices[:, :, 3]
    
    class_atkis_filtered = apply_multi_band_single_filter(raster, county_indices, atkis_indices)
    class_invekos_filtered = apply_multi_band_single_filter(raster, county_indices, invekos_indices)

    comb_filtered = cat(class_atkis_filtered, class_invekos_filtered; dims=Band)
    comb_filtered = setdims(comb_filtered, Band(1:2))
    return comb_filtered
end

```

and

```julia
function apply_multi_band_single_filter(prediction, county_index, inner_index)
    filtered_prediction = deepcopy(prediction)

    for u_county_idx in unique(county_index)
        county_mask = county_index .== u_county_idx
        coordinates = findall(county_mask)
    
        filter1d_inner_indices = inner_index[coordinates]
        Threads.@threads for u_inner_idx in unique(filter1d_inner_indices)
            tmp_inner_2d_mask = inner_index .== u_inner_idx
            tmp_inner_1d_mask = filter1d_inner_indices .== u_inner_idx

            temp_inner_mode = (u_inner_idx == 0) ? 0 : StatsBase.mode(filtered_prediction[coordinates][findall(tmp_inner_1d_mask)])

            filtered_prediction[findall(@. county_mask & tmp_inner_2d_mask)] .= temp_inner_mode
            end
        end

    return filtered_prediction
end

```

I split the indexes and calculated them seperately and than I wann to concatinate them again, to save both in one file.

The coordinates are from a common copied GeoArray, nevertheless writing the GeoArray to disc (as tif), results in this error:

> MethodError: no method matching indexorder(::AutoOrder)  
> Closest candidates are:  
> indexorder(!Matched::Tuple{}) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\dimension.jl:182  
> indexorder(!Matched::Ordered) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\mode.jl:38  
> indexorder(!Matched::Unordered) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\mode.jl:56  
> …
> 
> Stacktrace:  
> [1] indexorder(::Band{UnitRange{Int64},AutoMode{AutoOrder},Nothing}) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\dimension.jl:167  
> [2] reorderindex(::GeoArray{UInt8,3,Tuple{Lon{LinRange{Float64},Projected{Ordered{DimensionalData.Forward,DimensionalData.Forward,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Lat{LinRange{Float64},Projected{Ordered{DimensionalData.Reverse,DimensionalData.Reverse,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Band{UnitRange{Int64},AutoMode{AutoOrder},Nothing}},Tuple{},Array{UInt8,3},String,GDALarrayMetadata{String,Any},UInt16}, ::Band{DimensionalData.Forward,AutoMode{AutoOrder},Nothing}, ::DimensionalData.Forward) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\utils.jl:157  
> [3] reorderindex(::GeoArray{UInt8,3,Tuple{Lon{LinRange{Float64},Projected{Ordered{DimensionalData.Forward,DimensionalData.Forward,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Lat{LinRange{Float64},Projected{Ordered{DimensionalData.Reverse,DimensionalData.Reverse,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Band{UnitRange{Int64},AutoMode{AutoOrder},Nothing}},Tuple{},Array{UInt8,3},String,GDALarrayMetadata{String,Any},UInt16}, ::Band{DimensionalData.Forward,AutoMode{AutoOrder},Nothing}) at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\utils.jl:148  
> [4] reorderindex at C:\Users\TITAN.julia\packages\DimensionalData\IHaNW\src\utils.jl:144 [inlined]  
> [5] (::GeoData.var"#161#163")(::GeoArray{UInt8,3,Tuple{Lon{LinRange{Float64},Projected{Ordered{DimensionalData.Forward,DimensionalData.Forward,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Lat{LinRange{Float64},Projected{Ordered{DimensionalData.Reverse,DimensionalData.Reverse,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Band{UnitRange{Int64},AutoMode{AutoOrder},Nothing}},Tuple{},Array{UInt8,3},String,GDALarrayMetadata{String,Any},UInt16}) at C:\Users\TITAN.julia\packages\GeoData\MVllE\src\sources\gdal.jl:135  
> [6] |\> at .\operators.jl:834 [inlined]  
> [7] write(::String, ::Type{GDALarray}, ::GeoArray{UInt8,3,Tuple{Lon{LinRange{Float64},Projected{Ordered{DimensionalData.Forward,DimensionalData.Forward,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Lat{LinRange{Float64},Projected{Ordered{DimensionalData.Reverse,DimensionalData.Reverse,DimensionalData.Forward},Regular{Float64},Intervals{Start},WellKnownText{GeoFormatTypes.CRS,String},Nothing},GDALdimMetadata{Any,Any}},Band{UnitRange{Int64},AutoMode{AutoOrder},Nothing}},Tuple{},Array{UInt8,3},String,GDALarrayMetadata{String,Any},UInt16}) at C:\Users\TITAN.julia\packages\GeoData\MVllE\src\sources\gdal.jl:134  
> [8] top-level scope at In[36]:1  
> [9] include\_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 14, 2020, 11:54am UTC](https://discourse.julialang.org/t/using-archgdal/46006/11 "2020-09-14T11:54:42Z")

</div>

Ok this part is an easy fix

> [@Ordens\_Ritter](#):
>
> `setdims(comb_filtered, Band(1:2))`

You need to use:

```julia
setdims(comb_filtered, Band(1:2, NoIndex()))

```

To set the index type. Its annoying I know.  
The new `set` method fixes having to do any of that, but it isn’t finished yet.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 14, 2020, 12:06pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/12 "2020-09-14T12:06:57Z")

</div>

This isn’t the fastest way in julia:

```julia
raster[raster .== 55537] .= 0
raster[raster .== -9999] .= 0

```

You are allocating new geoarrays of booleans to index with. I’ll look into anything makes that especially slow with GeoArray compared to a regular array, but it’s better to do:

```julia
replace!(raster, 55537 => 0)
replace!(raster, -9999 => 0)

```

This just runs a function over every value in the array and updates it in place, so no allocations are made. Avoiding allocations is one of the main strategies for high performance julia.

But I guess the time is probably in those filter calls too.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 15, 2020, 1:38am UTC](https://discourse.julialang.org/t/using-archgdal/46006/13 "2020-09-15T01:38:53Z")

</div>

Things like this are allocating:

```julia
filtered_prediction[findall(@. county_mask & tmp_inner_2d_mask)] .= temp_inner_mode
filtered_prediction[coordinates][findall(tmp_inner_1d_mask)]

```

And your:

```julia
tmp_inner_1d_mask
tmp_inner_2d_mask

```

Could be defined outside the loop and reused (one for each thread) so you don’t keep allocating it.

---

<div class="post-metadata">

**Author:** ![Ordens\_Ritter](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ordens_ritter/32/50363_2.png) [@Ordens\_Ritter](https://discourse.julialang.org/u/Ordens_Ritter)\
**Post date:** [September 16, 2020, 3:10pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/14 "2020-09-16T15:10:22Z")

</div>

Thanks, @Raf

I think defining the type in the function with e. g.

local varname::GeoArray{2}

didn’t increase the speed, also it didn’t the accumulation. I am wondering, whether I could rearrange the code, that that I only wouldnneed a 2d or 1d array.

But because auf the speed issue I am reworking the algorithm. Therefore I want to ask, how to deal with the GeoArray, in more detail.

How do I apply a function like unique, maximum, find max for example on a GeoArray with a single Band or Multiple bands?

How can I extract the value of the bands, from the arrays Index, or how can I find the geo position from the index?

Perhaps it would be helpful,to add some usage examples for these on github,etc.  
Thanks a lot!

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [September 16, 2020, 5:43pm UTC](https://discourse.julialang.org/t/using-archgdal/46006/15 "2020-09-16T17:43:48Z")

</div>

If you want examples checkout the [docs for DimensionalData.jl](https://rafaqz.github.io/DimensionalData.jl/stable/). `GeoArray` is just a special case of `DimArray` from DimensionalData.jl.  
You can use maximum etc just like you would on a regular array, but you can use `Lat/Lon` etc as the `dims` argument if you need to.

I think your main problem for performance is the allocations you are making. `GeoArray` will be slower than normal array if you allocate or slice because it has the dim index as well. But in high performance code you shouldn’t do that anyway - and if you don’t slice or allocate GeoArray should be identical in performance to a regular array.

You should work out what is actually slow before you do anything else. Use ProfileView.jl to do this, or `@profiler` in Atom. And use BenchmarkTools.jl to time everything. Try to get you allocations to decrease and time every change.

Things like this don’t do anything - It will be type stable anyway:  
local varname::GeoArray{2}

Edit: also remember you can just get the parent Array from a GeoArray with `A = parent(geoarray)` and work with that, and copy your result back to a `GeoArray`, or use `geoarray = rebuild(geoarray; data=A)` to rebuild with new data (the same size of course).