# How to rasterize a geopackage (Corine Land Cover)?

**URL:** <https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316>\
**Category:** Geo\
**Created:** [June 7, 2024, 9:43am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316 "2024-06-07T09:43:26Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 9:43am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/1 "2024-06-07T09:43:26Z")

</div>

Hello, being only marginally on the domain, I am a bit lost with the various GIS packages…  
I need to rasterize a geopackage (gpkg) vector file (Corine Land Cover).  
Visualizing the gpkg file in QGis I can see there are various layers in a single file (in this case the main one is the European land use cover, the one I am interested in, and the others refer to overseas French territories).  
My final objective is to take the rasterization at high resolution and count the various pixels of a given class within a low resolution pixel, so I have as final product a set of rasters, each one indicating the percentual of land use for a specific class within each pixel.

I think I can do the second step, but I need first to do the first rasterization. I have tried in qgis to trsnform the geopackage layer as Shapefile, for which there is an example of Rasterization in `Rasters.jl`, but QGis complain that the file is too big for the Shapefile format.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 10:09am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/2 "2024-06-07T10:09:27Z")

</div>

An alternative approach that needs more testing is to load the geopackage directly with GeoIO.jl and then use the `Rasterize` transform from GeoStats.jl

If you encounter issues or performance problems, please report.

---

<div class="post-metadata">

**Author:** ![Fliks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fliks/32/2494_2.png) [@Fliks](https://discourse.julialang.org/u/Fliks)\
**Post date:** [June 7, 2024, 10:57am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/3 "2024-06-07T10:57:28Z")

</div>

You should be able to open the geopackage using ArchGDAL and then it will be the same as the shapefile example for rasterize. There is no need to convert the data beforehand, because both packages adhere to the GeoInterface and therefore these geo types are interchangeable.

---

<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:** [June 7, 2024, 11:13am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/4 "2024-06-07T11:13:42Z")

</div>

Yes ArchGDAL or GeoDataframes.jl will open that fine, or also GeoIO.jl.

Then Rasters.jl should just rasterize any of those objects directly in `rasterize` without modification, because GeoInterface.jl handles the conversions for you.

There are no comprehensive benchmarks at this stage but it looks like `Rasters.rasterize` its one of the fastest implementations that exists. If you find anything else as fast, make an issue 😉

(Also we started on a native julia [GitHub - JuliaGeo/GeoPackage.jl: A Julia reader and writer for `.gpkg` files](https://github.com/JuliaGeo/GeoPackage.jl) but didn’t get so far yet)

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 11:40am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/5 "2024-06-07T11:40:19Z")

</div>

Thank you both @Fliks and @Raf for the quick answers.

So, trying to work on your answers. The file I am trying to read is the Corine Land Cover 2018, in geopackage format.

- official: [https://land.copernicus.eu/en/products/corine-land-cover/clc2018](https://land.copernicus.eu/en/products/corine-land-cover/clc2018)
- doi: [https://doi.org/10.2909/71c95a07-e296-44fc-b22b-415f42acfdf0](https://doi.org/10.2909/71c95a07-e296-44fc-b22b-415f42acfdf0)
- downlodable copy on our server: [https://nc.beta-lorraine.fr/s/9NrqQBmR5ib2wmd/download](https://nc.beta-lorraine.fr/s/9NrqQBmR5ib2wmd/download)

```julia
using Downloads, Rasters, GeoIO, ArcGDAL
url = "https://nc.beta-lorraine.fr/s/9NrqQBmR5ib2wmd/download"
clc2018_path = "clc2018.gpkg"
Downloads.download(url,clc2018_path)
table = GeoIO.load(clc2018_path) # ERROR: geometry column not found
table = ArchGDAL.read(clc2018_path)

```

The ArchGDAL call produces the following object (superquickly, so it is a lazy object):

```julia
julia> table
GDAL Dataset (Driver: GPKG/GeoPackage)
File(s): 
  /home/lobianco/CloudFiles/beta-lorraine-sync/Documents/DatiGeografici/EU/corine/clc2018/u2018_clc2018_v2020_20u1_geoPackage/DATA/U2018_CLC2018_V2020_20u1.gpkg

Number of feature layers: 6
  Layer 0: U2018_CLC2018_V2020_20u1 (wkbMultiPolygon)
  Layer 1: U2018_CLC2018_V2020_20u1_FR_REU (wkbMultiPolygon)
  Layer 2: U2018_CLC2018_V2020_20u1_FR_GLP (wkbMultiPolygon)
  Layer 3: U2018_CLC2018_V2020_20u1_FR_GUF (wkbMultiPolygon)
  Layer 4: U2018_CLC2018_V2020_20u1_FR_MTQ (wkbMultiPolygon)
  Remaining layers:
    U2018_CLC2018_V2020_20u1_FR_MYT, 

```

However, I then have an error in Rasters.raster():

```julia
clc_raster = Rasters.rasterize(last, table,fill=1) # ERROR: MethodError: no method matching iterate(::ArchGDAL.IDataset)
clc_raster = Rasters.rasterize(last, table[1],fill=1) # ERROR: MethodError: no method matching getindex(::ArchGDAL.IDataset, ::Int64)

```

Any clue ?

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 11:54am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/6 "2024-06-07T11:54:44Z")

</div>

I managed to get the layer, but still it seems the ArchGDAL object is not compatible:

```julia
layer1 = ArchGDAL.getlayer(table, 0)
clc_raster = Rasters.rasterize(last, layer1,res=10000, missingval=0, fill=1, progress=true) 

```

Results in :

`ERROR: ArgumentError: Object is not a GeoInterface.jl compatible geometry: ()`

---

<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:** [June 7, 2024, 11:56am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/7 "2024-06-07T11:56:16Z")

</div>

Hm seems these objects are not identifying as a GeoInterface.jl geometry/feature/feature collection or with a Tables.jl table.

Rasters.jl doesn’t actually know about ArchGDAL objects, it just follows the Tables.jl and GeoInterface.jl interfaces.

Try using GeoDataFrames.jl instead? That should handle it for you. ArchGDAL is pretty low level.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 11:56am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/8 "2024-06-07T11:56:43Z")

</div>

Thanks @juliohm but `GeoIO.load(filepath)` results in an `ERROR: geometry column not found` ☹

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 11:57am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/9 "2024-06-07T11:57:47Z")

</div>

Thank you for reporting the error @sylvaticus , it is probably in the file itself. We will investigate and report our findings.

Notice that, at present, GeoIO.jl simply relies on ArchGDAL.jl for .gpkg:

> <https://github.com/JuliaEarth/GeoIO.jl/blob/7865d065b1d3b9f4727d99b97988dc0bfa874840/src/load.jl#L83-L94>

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 11:59am UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/10 "2024-06-07T11:59:45Z")

</div>

thanks… to note, I can open it with qgis without issues…

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 12:02pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/11 "2024-06-07T12:02:29Z")

</div>

I’ve updated my comment above with more info so that you can inspect the issue in more depth, in parallel.

My guess is that the file itself (still downloading here) lacks the expected structure for the backend packages to load it.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 12:09pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/12 "2024-06-07T12:09:58Z")

</div>

Ok, apparently the `geometry` column is named `Shape` in this file, which is not covered yet:

> <https://github.com/JuliaEarth/GeoIO.jl/blob/7865d065b1d3b9f4727d99b97988dc0bfa874840/src/utils.jl#L18-L28>

Let me try to work on a hot fix for this.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 12:12pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/13 "2024-06-07T12:12:44Z")

</div>

yes, it is a multi-vector layer.  
The `read` call of ArchGDAL just get the whole package, but then one has to use, as I discovered, use `ArchGDAL.getlayer()` to get the specific layer.

I will try 2 more attempts now:

1. use GeoDataframes.jl to load the layer
2. trying to use the Rasterize transform from GeoStats.jl on the `ArchGDAL.getlayer` output…

Edit:  
First attempt succeeded in creating a df (got a warning " Warning: This file has multiple layers, you only get the first layer by default now." but that’s fine for me).  
I obtained a df that looks like:

```julia
julia> table
2375406×5 DataFrame
     Row │ Shape Code_18 Remark Area_Ha ID         
         │ IGeometr… String String? Float64 String     
─────────┼───────────────────────────────────────────────────────────────────
       1 │ Geometry: wkbMultiPolygon 111 missing 130.864 EU_1
       2 │ Geometry: wkbMultiPolygon 111 missing 53.7445 EU_2
       3 │ Geometry: wkbMultiPolygon 111 missing 30.7191 EU_3
       4 │ Geometry: wkbMultiPolygon 111 missing 50.2018 EU_4
       5 │ Geometry: wkbMultiPolygon 111 missing 481.849 EU_5
       6 │ Geometry: wkbMultiPolygon 111 missing 42.6509 EU_6
    ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮
 2375402 │ Geometry: wkbMultiPolygon 512 missing 807.041 EU_2375402
 2375403 │ Geometry: wkbMultiPolygon 512 missing 141.363 EU_2375403
 2375404 │ Geometry: wkbMultiPolygon 512 missing 246.634 EU_2375404
 2375405 │ Geometry: wkbMultiPolygon 512 missing 41.8458 EU_2375405
 2375406 │ Geometry: wkbMultiPolygon 512 missing 43.3464 EU_2375406
                                                         2375395 rows omitted

```

However passing this df to `Rasters.rasterize` I got the same error as the gdal object (`ERROR: ArgumentError: Object is not a GeoInterface.jl compatible geometry: ()`)

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 12:27pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/14 "2024-06-07T12:27:07Z")

</div>

Ok, prepared a hotfix to handle this unexpected geometry column name:

```julia
# helper function to find the
# geometry column of a table
function geomcolumn(names)
  snames = string.(names)
  gnames = ["geometry","geom","shape"]
  gnames = [gnames; uppercasefirst.(gnames)]
  gnames = [gnames; uppercase.(gnames)]
  gnames = [gnames; [""]]
  select = findfirst(∈(snames), gnames)
  if isnothing(select)
    throw(ErrorException("geometry column not found"))
  else
    Symbol(gnames[select])
  end
end

```

Managed to load the file without problems with the master branch of GeoIO.jl:

```julia
julia> GeoIO.load("/home/juliohm/clc2018.gpkg")
2375406×5 GeoTable over 2375406 GeometrySet
┌─────────────┬─────────────┬────────────┬─────────────┬───────────────────┐
│ Code_18 │ Remark │ Area_Ha │ ID │ geometry │
│ Categorical │ Categorical │ Continuous │ Categorical │ MultiPolygon │
│ [NoUnits] │ [NoUnits] │ [NoUnits] │ [NoUnits] │ │
├─────────────┼─────────────┼────────────┼─────────────┼───────────────────┤
│ 111 │ missing │ 130.864 │ EU_1 │ Multi(1×PolyArea) │
│ 111 │ missing │ 53.7445 │ EU_2 │ Multi(1×PolyArea) │
│ 111 │ missing │ 30.7191 │ EU_3 │ Multi(1×PolyArea) │
│ 111 │ missing │ 50.2018 │ EU_4 │ Multi(1×PolyArea) │
│ 111 │ missing │ 481.849 │ EU_5 │ Multi(1×PolyArea) │
│ 111 │ missing │ 42.6509 │ EU_6 │ Multi(1×PolyArea) │
│ 111 │ missing │ 139.015 │ EU_7 │ Multi(1×PolyArea) │
│ 111 │ missing │ 76.7515 │ EU_8 │ Multi(1×PolyArea) │
│ 111 │ missing │ 308.071 │ EU_9 │ Multi(1×PolyArea) │
│ 111 │ missing │ 127.451 │ EU_10 │ Multi(1×PolyArea) │
│ 111 │ missing │ 84.7652 │ EU_11 │ Multi(1×PolyArea) │
│ ⋮ │ ⋮ │ ⋮ │ ⋮ │ ⋮ │
└─────────────┴─────────────┴────────────┴─────────────┴───────────────────┘
                                                        2375395 rows omitted

```

You can forward the `layer` option if necessary with `GeoIO.load(file, layer=1)`.

Triggered a patch release so that you can update your environment:

> <https://github.com/JuliaRegistries/General/pull/108455>
>
> \- Registering package: GeoIO
> \- Repository: https://github.com/JuliaEarth/GeoIO.j…l
> \- Created by: @juliohm
> \- Version: v1.13.1
> \- Commit: 2a4e1cac6003a23e4004e808835f78f7501206e8
> \- Reviewed by: @juliohm
> \- Reference: https://github.com/JuliaEarth/GeoIO.jl/commit/2a4e1cac6003a23e4004e808835f78f7501206e8#commitcomment-142846458
> \- Description: Load/save geospatial data compatible with the GeoStats.jl framework

---

<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:** [June 7, 2024, 12:28pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/15 "2024-06-07T12:28:23Z")

</div>

Ok one problem is the geometry column name is weird. You can manually pass that to Rasters.jl using the keyword `geometrycolumn=:Shape`

Or just get that vector out of the data frame manually and rasterize it.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 12:31pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/16 "2024-06-07T12:31:24Z")

</div>

@sylvaticus if you decide to try the `Rasterize` transform on the loaded geotable with GeoIO.jl, please feel free to report further issues. These reports really help improve robustness of the packages.

---

<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:** [June 7, 2024, 12:34pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/17 "2024-06-07T12:34:24Z")

</div>

If people want to use another package that’s fine, but honestly I find these threads with two competing methods from different packages are pretty distracting.

The question was for Rasters.jl and the solution is very simple. Let’s not do this all the time. We can wait to see if there is no easy solution to the direct question.

In this case there is a keyword for exactly this problem.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [June 7, 2024, 12:36pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/18 "2024-06-07T12:36:43Z")

</div>

I wouldn’t join the thread with an alternative solution if the title was requesting Rasters.jl specifically. I thought that @sylvaticus just wanted to accomplish a generic task as a first-time user of GIS software.

If this is a Rasters.jl question, sorry for the confusion.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [June 7, 2024, 12:44pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/19 "2024-06-07T12:44:20Z")

</div>

I am very sorry, I didn’t want to spark an issue between GIS packages.  
Just GIS is not my main domain, but I needed to do “some GIS” within a larger model.  
I have manually renamed the df outputed by GeoDataFrames and this then indeed worked with Rasters.Rasterize(), I will also try the GeoIO → GeoStats.Rasterize method and report if it works…

I have a related issue but I am now afraid to ask… :-/ How do I rasterize by setting the pixel value the one with the majority land use in the pixel (land use is in the vector format in the `Code_18` column).

---

<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:** [June 7, 2024, 12:48pm UTC](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316/20 "2024-06-07T12:48:39Z")

</div>

Op mentioned Rasters.jl in main comment and then shows a broken Rasters.rasterize example in the next comment… I’m just here trying to make that work.

@sylvaticus no worries at all

What do you mean the majority land use? Do you mean here multiple geometries cover the same pixel?

[Next page](https://discourse.julialang.org/t/how-to-rasterize-a-geopackage-corine-land-cover/115316.md?page=2)