# Retrieving CRS from vector and raster data

**URL:** <https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903>\
**Category:** Geo\
**Tags:** geometry, polygons, rasters, geodesy\
**Created:** [December 21, 2023, 12:55pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903 "2023-12-21T12:55:53Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Bojas](https://avatars.discourse-cdn.com/v4/letter/b/8e7dd6/32.png) [@Bojas](https://discourse.julialang.org/u/Bojas)\
**Post date:** [December 21, 2023, 12:55pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/1 "2023-12-21T12:55:53Z")

</div>

Hello, I’m new to Julia, as for now I was using mainly Python. What I want to do is to crop part of raster based on polygon. Both polygon and raster are in different spatial references and I have to convert either one of this (preferably polygon). The problem is I don’t know how I can retrieve CRS from polygon and raster objects. In python I used rasterio library to read raster CRS and pyproj library to convert polygon geometry held in GeoDataFrame. In Julia, I see there is a PROJ module, but I cannot find a way to retrieve raster’s CRS, even if I see crs param in output.

```julia
import GeoDataFrames as GDF
using Rasters

gdf = GDF.read("./geometry/warsaw_3857.shp") # GeoDataFrame containing one polygon; As stated in the name, it's polygon of Warsaw in EPSG:3857

geom = gdf.geometry[1]

my_raster = Raster("./rasters/red.jp2") # raster with red band from Sentinel-2 in EPSG:32634

```

Is there a simple way to obtain raster’s CRS in Julia? Should I use another packages to read data?

---

<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:** [December 21, 2023, 12:59pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/2 "2023-12-21T12:59:03Z")

</div>

`Rasters.crs(raster)` 😉

for the `gdf` you may need `GDF.crs(gdf)` ? or `GeoInterface.crs` ?

This is actually a good reminder to remove `Rasters.crs` and instead use `GeoInteface.crs` - I just never got around to it.

---

<div class="post-metadata">

**Author:** ![Bojas](https://avatars.discourse-cdn.com/v4/letter/b/8e7dd6/32.png) [@Bojas](https://discourse.julialang.org/u/Bojas)\
**Post date:** [December 21, 2023, 1:17pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/3 "2023-12-21T13:17:05Z")

</div>

Thank you, Rasters.crs(my\_raster) worked successfuly, although it returned long xml-like string:

> WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), “PROJCS["WGS 84 [/](https://vscode-remote+ssh-002dremote-002b64-002e225-002e139-002e4.vscode-resource.vscode-cdn.net/) UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS\_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse\_Mercator"],PARAMETER["latitude\_of\_origin",0],PARAMETER["central\_meridian",21],PARAMETER["scale\_factor",0.9996],PARAMETER["false\_easting",500000],PARAMETER["false\_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32634"]]”)

unlike rasterio in python, which returned exact crs string - “EPSG:32634”, but I guess it can’t be helped. What I would like to ask is why GeoInterface.crs(geom) does not return anything. In shapefile CRS is declared so I wonder why Julia can’t see it? Besides nor GDF.crs(gdf) nor GeoInterface.crs(my\_raster) seems to work, or I’m doing something wrong. After using either one of the stated method I received this error:

> MethodError: no method matching crs(::Nothing, ::Raster{UInt16, 2, Tuple{X{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, X{Colon}}}, Y{Projected{Float64, LinRange{Float64, Int64}, DimensionalData.Dimensions.LookupArrays.ReverseOrdered, DimensionalData.Dimensions.LookupArrays.Regular{Float64}, DimensionalData.Dimensions.LookupArrays.Intervals{DimensionalData.Dimensions.LookupArrays.Start}, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, WellKnownText{GeoFormatTypes.CRS}, Nothing, Y{Colon}}}}, Tuple{Band{DimensionalData.Dimensions.LookupArrays.Categorical{Int64, UnitRange{Int64}, DimensionalData.Dimensions.LookupArrays.ForwardOrdered, DimensionalData.Dimensions.LookupArrays.NoMetadata}}}, Matrix{UInt16}, Symbol, DimensionalData.Dimensions.LookupArrays.Metadata{Rasters.GDALsource, Dict{String, Any}}, Nothing})

---

<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:** [December 21, 2023, 2:27pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/4 "2023-12-21T14:27:00Z")

</div>

That “XML” is Well Known Text, the default crs standard in GDAL.

Conversion of Well Known Text to EPSG numbers is not 100% accurate or always possible… so I’m not sure why rasterio would do that by default. I guess it does look nicer.

`GeoInterface.crs` does not work on a Raster yet, as I mentioned I have forgotten to update it as the Rasters method existed first. I’ll change that soon and make them use the same method (this quesiton is a good reminder to do that)

Currently `GeoInterface.crs` will work on shapefile/polygons etc.

This PR will make them the same method:

> <https://github.com/rafaqz/Rasters.jl/pull/580>

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [December 21, 2023, 3:38pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/5 "2023-12-21T15:38:49Z")

</div>

> [@Bojas](#):
>
> ```julia
> gdf = GDF.read("./geometry/warsaw_3857.shp") # GeoDataFrame containing one polygon; As stated in the name, it's polygon of Warsaw in EPSG:3857
> 
> geom = gdf.geometry[1]
> 
> ```

Ok, loads of bug reports to make, but here’s how it should work and what works now. Edit: @Bojas Releases have been made, if you update your packages this should work:

```julia
import GeoFormatTypes as GFT
using GeoInterface

GeoInterface.crs(geom)

```

However, the crs is also stored on the DataFrame level, in the metadata.

```julia
GDF.metadata(df)["crs"]

```

---

<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:** [December 21, 2023, 4:21pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/6 "2023-12-21T16:21:42Z")

</div>

What do you think `GeoInterface.crs(geom)` should return?  
Should this simply be a call of `AG.getspatialref(geom)`?

Should the conversion of the AG.getspatialref return tothe GeoInterfaceTypes type be necessary or could we make the AG type a subtype of GeoInterfaceTypes?

---

<div class="post-metadata">

**Author:** ![evetion](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/evetion/32/22679_2.png) [@evetion](https://discourse.julialang.org/u/evetion)\
**Post date:** [December 21, 2023, 4:35pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/7 "2023-12-21T16:35:15Z")

</div>

> What do you think `GeoInterface.crs(geom)` should return?  
> Should this simply be a call of `AG.getspatialref(geom)`?

> Should the conversion of the AG.getspatialref return tothe GeoInterfaceTypes type be necessary or could we make the AG type a subtype of GeoInterfaceTypes?

As document by GeoInterface, a GeoFormatTypes thing. And yes, that should be a call to getspatialref, currently making a PR.

Not sure what you mean by GeoInterfaceTypes, if its GeoFormatTypes, I don’t think the AbstractSpatialRefs would be a good subset of GeoFormatTypes. Although we could do with an `iscrs` trait?

---

<div class="post-metadata">

**Author:** ![JosephPollacco](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josephpollacco/32/51404_2.png) [@JosephPollacco](https://discourse.julialang.org/u/JosephPollacco)\
**Post date:** [April 1, 2025, 2:15pm UTC](https://discourse.julialang.org/t/retrieving-crs-from-vector-and-raster-data/107903/8 "2025-04-01T14:15:00Z")

</div>

> [@evetion](#):
>
> AG.getspatialref(geom

This cam be an alternative solution:

```julia
Grid_GeoTIFF = GeoTIFF.load(Path)
            Grid_GeoTIFF_Metadata = GeoTIFF.metadata(Grid_GeoTIFF)
            Crs =GeoTIFF.epsgcode(Grid_GeoTIFF_Metadata) |>Int

```