GeoTiff files retrieving metadata by using Rasters.jl as in rasterio.py

Specific Domains Geo
, ,
1

Hi,

I am looking how to retrive metadata from GeoTiff files by using Julia language.

I can obtain the relevant information by using the rasterio.jl python package, but I am struggling to perform the same commands by using the Julia language. I looked at the Rasters.jl package, but could not find an easy answer.

Below is the code I wrote in Python and it will be great to find similar package in Julia language.

Many thanks for any help you may provide,
Joseph

import rasterio

def INFORMATION(Path_Input, Verbose=True):

	Rasterio = rasterio.open(Path_Input)
	Dimensions   = Rasterio.transform

	CellSize_X   = Dimensions[0]
	CellSize_Y   = Dimensions[4]
	N_Width      = Rasterio.width
	N_Hight      = Rasterio.height
	CRS          = Rasterio.crs
	Coord_X_Left   = Rasterio.bounds[0]
	Coord_X_Right  = Rasterio.bounds[2]
	Coord_Y_Top    = Rasterio.bounds[3]
	Coord_Y_Bottom = Rasterio.bounds[1]

	Coord_TopLeft     = [Coord_Y_Top, Coord_X_Left]
	Coord_BottomRight = [Coord_Y_Bottom, Coord_X_Right]

	if Verbose:
		print("Cell size_X=", CellSize_X )
		print("Cell size_Y=", CellSize_Y )
		print("DEM width= " , N_Width) 
		print("DEM height= ", N_Hight)
		print("Coordinate Reference System CRS=" , CRS)
		print("Bounds= ", Rasterio.bounds)
		print("Bound_Left= ", Coord_X_Left ," ,Bound_Right= ", Coord_X_Right ," ,Bound_Top= ", Coord_Y_Top ," ,Bound_Botton= ", Coord_Y_Bottom)
		print("Coord_TopLeft= ", Coord_TopLeft)
		print("Coord_BottomRight= ", Coord_BottomRight)

	return CellSize_X, CellSize_Y, Coord_BottomRight, Coord_TopLeft, Coord_X_Left, Coord_X_Right, Coord_Y_Bottom, Coord_Y_Top, CRS, N_Hight, N_Width, Rasterio
2

With GMT.jl you’do (e.g.)

meta = gdalinfo("C:\\v\\L9\\RGB_surf.tiff");

println(meta)
Driver: GTiff/GeoTIFF
Files: C:\v\L9\RGB_surf.tiff
Size is 7731, 7841
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 29N",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            MEMBER["World Geodetic System 1984 (G2296)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 29N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Navigation and medium accuracy spatial referencing."],
        AREA["Between 12Β°W and 6Β°W, northern hemisphere between equator and 84Β°N, onshore and offshore. Algeria. CΓ΄te D'Ivoire (Ivory Coast). Faroe Islands. Guinea. Ireland. Jan Mayen. Liberia, Mali. Mauritania. Morocco. Portugal. Sierra Leone. Spain. United Kingdom (UK). Western Sahara."],
        BBOX[0,-12.01,84.01,-6]],
    ID["EPSG",32629]]
Data axis to CRS axis mapping: 1,2
Origin = (394185.000000000000000,4423215.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  394185.000, 4423215.000) ( 10d14'19.54"W, 39d57' 8.84"N)
Lower Left  (  394185.000, 4187985.000) ( 10d12' 8.92"W, 37d49'59.37"N)
Upper Right (  626115.000, 4423215.000) (  7d31'25.17"W, 39d56'58.83"N)
Lower Right (  626115.000, 4187985.000) (  7d34' 0.82"W, 37d49'50.08"N)
Center      (  510150.000, 4305600.000) (  8d52'58.62"W, 38d53'57.21"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
1 Like
3

You can use Rasters to derive the information you are looking for:

r = Raster(somepath, lazy=true)
CellSize_X = step(dims(r, X))
Cellsize_Y = step(dims(r, Y))
N_Width = size(r, X)
N_heigth = size(r, Y)
CRS = crs(r)
Coord_X_Left = first(dims(r,X))
Coord_X_Right = last(dims(r,X))
Coord_Y_Top = first(dims(r,Y))
Coord_Y_Bottom = last(dims(r,Y))

It might be, that the Coord values are switched depending on the orientation of the dimension. You would have to lookout for ReverseOrdered in the printout of the dimension.

Where did you look for the info and is there anything that we could improve our docs on so that you would have found the information easier?

4

Thanks @joa-quim that is so kind of you. I will need to automatically check the metadata. Now the question would be how can you automatically put into parameters for e.g. Upper Left, Lower Left or Pixel Size etc…?

5

@JosephPollacco the native GeoTIFF.jl reader/writer has all the metadata you need.

If you just need to load the data correctly with it, look into GeoIO.jl.

The packages listed above rely on the external GDAL library. You can decide which option is more adequate for your use case.

6

That is really amazing and a great thanks for your response. I tested it and it all works amazing.

A quick question with CRS how can I query the EPSG value?

For the improvment of the already excellent Rasters documentation , it might be great that you include the following example for retrieving metadata, and I am wandering if you could make simple functions following Rasterio e.g.

Raster.CellSize_X(…) , Raster.Extent(…)

Many thanks for your great work,
Joseph

7

@juliohm Thanks for your usefull information, not sure how to abstract the metadata of interest with GeoTIFF.jl

8

@JosephPollacco if you are familiar with the GeoTIFF specification, you will be able to access any metadata stored in the file. First you retrieve the metadata object from the loaded image:

image = GeoTIFF.load(filename)

metadata = GeoTIFF.metadata(image)

and then you query the keys of interest:

# query the raster type
GeoTIFF.geokey(metadata, GTRasterTypeGeoKey)

You can find the list of all keys here:

If the information you need is not directly available as a key, but instead is derived from the keys above, it shouldn’t be hard to write a small helper function to get the value you need.

Keep in mind that GeoTIFF.jl adheres to the GeoTIFF specification strictly. The GDAL reader does more things with GDAL-specific keys (not part of the GeoTIFF specification). Try to not rely on these GDAL-specific keys if you can.

9

Actually we have Raaters.extent(raster) for the extent.

We don’t always know cellsize as NetCSDF etc may have arbitrary/mixed sizes along a dimension. But maybe we should have a function that can return nothing for some axes.

But you can do map(step, dims(raster)), and it will just fail for Irregular dimension as they have no step size.

10

With risk of this getting very confusing with the abundance of alternatives, and quickly as I have to leave for the whole afternoon. Other options with GMT would be (ignore the z_min/z_max that are garage when the array is of UInt8)

julia> grdinfo("C:\\v\\L9\\RGB_surf.tiff", C=true)

1Γ—12 GMTdataset{Float64, 2}
 Row β”‚    x_min     x_max      y_min      y_max        z_min         z_max    dx    dy  n_cols  n_rows  reg  isgeog
─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 β”‚ 394185.0  626115.0  4.18798e6  4.42322e6  1.79769e308  -1.79769e308  30.0  30.0  7731.0  7841.0  1.0     0.0

GMT.getproj("C:\\v\\L9\\RGB_surf.tiff")
"PROJCS[\"WGS 84 / UTM zone 29N\",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\",-9],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\",\"32629\"]]"

or

julia> GMT.getproj("C:\\v\\L9\\RGB_surf.tiff", proj4=true)
"+proj=utm +zone=29 +datum=WGS84 +units=m +no_defs"
11

Also @JosephPollacco if you want any changes to Rasters.jl github issues are the way! We can totally add a function to make getting the step size easier, but without a github issue your problem doesn’t really exist in an operational sense - thats where triage on needed work happens whenever we have time to do it. So most of the time its best to go straight to a github issue with literally any problems at all.

And discourse discussions about Rasters.jl tend to generate a lot of unrelated noise :wink:

1 Like
12

We will do, thanks

13

Thanks, it is great that there are different packages for which we can derive the metadata. Based on the information, I am proposing a solution, please feel free to improve it:

using ArchGDAL
const AG = ArchGDAL
using Rasters, GeoTIFF, Base

    Base.@kwdef mutable struct METADATA
        N_Width        :: Int64
        N_Height       :: Int64
        Ξ”X             :: Int64
        Ξ”Y             :: Int64
        Coord_X_Left   :: Float64
        Coord_X_Right  :: Float64
        Coord_Y_Top    :: Float64
        Coord_Y_Bottom :: Float64
        Crs            :: Int64
        Bands          :: Int64
    end # struct METADATA
    """
        Deriving metadata from the GeoTiff file
    """
    # ================================================================
    #		FUNCTION : RASTER_METADATA
    # ================================================================
        function RASTER_METADATA(Path)

            Grid = Rasters.Raster(Path, lazy=true)
                N_Width = size(Grid, X)
                N_Height = size(Grid, Y)
                Ξ”X =  step(dims(Grid, X)) |> abs
                Ξ”Y =  step(dims(Grid, Y)) |> abs
                Coord_X_Left   = first(dims(Grid, X))
                Coord_X_Right  = last(dims(Grid, X))
                Coord_Y_Top    = first(dims(Grid ,Y))
                Coord_Y_Bottom = last(dims(Grid,Y))

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

            Grid_Ag = AG.readraster(Path)
                Bands = AG.nraster(Grid_Ag)

            Metadata = METADATA(N_Width, N_Height, Ξ”X, Ξ”Y, Coord_X_Left, Coord_X_Right,Coord_Y_Top, Coord_Y_Bottom, Crs, Bands)

        return Metadata
        end # function RASTER_METADATA 
    # ----------------------------------------------------------------