From/to Rasterio

Note

Go to the end to download the full example code.

From/to Rasterio#

This example demonstrates importing and exporting a rasterio.io.DatasetReader or rasterio.io.MemoryFile from and to an xarray.DataArray with the rst accessor.

import rasterio as rio

A raster can be imported from a rasterio.io.DatasetReader with open_raster().

import geoutils as gu

filename_rast = gu.examples.get_path("exploradores_aster_dem")
ds = rio.open(filename_rast)
rast = gu.open_raster(ds)
rast

<xarray.DataArray (y: 618, x: 539)> Size: 1MB
[333102 values with dtype=float32]
Coordinates:
  * y            (y) float64 5kB 4.852e+06 4.852e+06 ... 4.834e+06 4.834e+06
  * x            (x) float64 4kB 6.272e+05 6.272e+05 ... 6.433e+05 6.433e+05
    band         int64 8B 1
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:             Area
    scale_factor:              1.0
    add_offset:                0.0
    overview_factors:          []
    overview_selected_level:   -1
    overview_selected_factor:  1
    block_shape:               (3, 539)


The data is unloaded, as when opened with a filename. The data will be loaded explicitly by any function requiring its values, such as plotting.

rast.rst.plot(cmap="terrain")
import raster

We can also open a rasterio.io.MemoryFile.

mem = rio.MemoryFile(open(filename_rast, "rb"))
rast = gu.open_raster(mem)
rast

<xarray.DataArray (y: 618, x: 539)> Size: 1MB
[333102 values with dtype=float32]
Coordinates:
  * y            (y) float64 5kB 4.852e+06 4.852e+06 ... 4.834e+06 4.834e+06
  * x            (x) float64 4kB 6.272e+05 6.272e+05 ... 6.433e+05 6.433e+05
    band         int64 8B 1
    spatial_ref  int64 8B 0
Attributes:
    AREA_OR_POINT:             Area
    scale_factor:              1.0
    add_offset:                0.0
    overview_factors:          []
    overview_selected_level:   -1
    overview_selected_factor:  1
    block_shape:               (3, 539)


The raster has the same DataArray representation.

Finally, we can export the xarray.DataArray to a rasterio.io.MemoryFile and open it as a rasterio.io.DatasetReader.

with rio.MemoryFile() as output_memory:
    rast.rst.to_file(output_memory.name)
    with output_memory.open() as output_dataset:
        print(output_dataset.profile)

{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.0, 'width': 539, 'height': 618, 'count': 1, 'crs': CRS.from_wkt('PROJCS["WGS 84 / UTM zone 18S",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",-75],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","32718"]]'), 'transform': Affine(30.0, 0.0, 627175.0,
       0.0, -30.0, 4852085.0), 'blockxsize': 539, 'blockysize': 3, 'tiled': False, 'interleave': 'band'}

Total running time of the script: (0 minutes 0.541 seconds)

Gallery generated by Sphinx-Gallery