Retrieving CRS from vector and raster data

Specific Domains Geo
, , ,
1

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.

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?

1 Like
2

Rasters.crs(raster) :wink:

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.

3

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

WellKnownText{GeoFormatTypes.CRS}(GeoFormatTypes.CRS(), “PROJCS["WGS 84 / 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})

4

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:

1 Like
5

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:

import GeoFormatTypes as GFT
using GeoInterface

GeoInterface.crs(geom)

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

GDF.metadata(df)["crs"]
6

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?

7

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?