Reading and plotting GeoTiff with Xarray

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import matplotlib.colors as mcolors

dset = xr.open_dataset("srtm_41_19.tif", engine="rasterio")

dset

<xarray.Dataset>
Dimensions:      (band: 1, x: 837, y: 661)
Coordinates:
  * band         (band) int64 1
  * x            (x) float64 9.695e+05 9.696e+05 ... 1.038e+06 1.038e+06
  * y            (y) float64 6.25e+06 6.25e+06 6.25e+06 ... 6.196e+06 6.196e+06
    spatial_ref  int64 ...
Data variables:
    band_data    (band, y, x) float32 ...

xarray.Dataset

• Dimensions:
  • band: 1
  • x: 837
  • y: 661
• Coordinates: (4)
  • band
    (band)
    int64
    1

    array([1])

  • x
    (x)
    float64
    9.695e+05 9.696e+05 ... 1.038e+06

    array([ 969532.272117,  969614.265543,  969696.258969, ..., 1037914.789564,
           1037996.78299 , 1038078.776417])

  • y
    (y)
    float64
    6.25e+06 6.25e+06 ... 6.196e+06

    array([6250255.998851, 6250174.005425, 6250092.011999, ..., 6196304.324414,
           6196222.330988, 6196140.337562])

  • spatial_ref
    ()
    int64
    ...

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 33S",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",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32733"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    WGS 84 / UTM zone 33S

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    15.0

    false_easting :

    500000.0

    false_northing :

    10000000.0

    scale_factor_at_central_meridian :

    0.9996

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 33S",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",15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32733"]]

    GeoTransform :

    969491.2754036566 81.99342619588413 0.0 6250296.995564442 0.0 -81.99342619588413

    [1 values with dtype=int64]

• Data variables: (1)
  • band_data
    (band, y, x)
    float32
    ...

    AREA_OR_POINT :

    Area

    [553257 values with dtype=float32]

• Indexes: (3)
  • band
    PandasIndex

    PandasIndex(Index([1], dtype='int64', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([ 969532.2721167546,  969614.2655429505,  969696.2589691463,
            969778.2523953422,   969860.245821538,   969942.239247734,
            970024.2326739299,  970106.2261001257,  970188.2195263216,
            970270.2129525176,
           ...
           1037340.8355807507, 1037422.8290069466, 1037504.8224331425,
           1037586.8158593384, 1037668.8092855343, 1037750.8027117301,
            1037832.796137926, 1037914.7895641219, 1037996.7829903178,
           1038078.7764165137],
          dtype='float64', name='x', length=837))

  • y
    PandasIndex

    PandasIndex(Index([ 6250255.998851344,  6250174.005425148,  6250092.011998952,
            6250010.018572756,  6249928.025146561, 6249846.0317203645,
            6249764.038294169,  6249682.044867973,  6249600.051441777,
            6249518.058015581,
           ...
            6196878.278397824,  6196796.284971627,  6196714.291545432,
            6196632.298119236, 6196550.3046930395,  6196468.311266844,
            6196386.317840648,  6196304.324414453,  6196222.330988256,
            6196140.337562061],
          dtype='float64', name='y', length=661))

• Attributes: (0)

Visualize (default QGIS)

dset["band_data"].plot(cmap="Greys_r")

<matplotlib.collections.QuadMesh at 0x7f1f16358250>

Visualize with pseudo colour

rgb_colors = np.array([[  0.,  80., 191.],
 [ 42., 152. ,  8.],
 [255. ,248. , 57.],
 [255. ,166. ,  0.],
 [189. , 11.,  35.]])

rgb_colors = rgb_colors/255.
my_colormap = mcolors.LinearSegmentedColormap.from_list('colormap', rgb_colors)

dset["band_data"].plot(cmap=my_colormap)

<matplotlib.collections.QuadMesh at 0x7f1f163d3bd0>