Opening GeoTIFF and NetCDF files with xarray ¶

**Sign up to the DEA Sandbox** to run this notebook interactively from a browser
Compatibility: Notebook currently compatible with both the NCI and DEA Sandbox environments

Background¶

It can be useful open an external raster dataset that you have previously saved to GeoTIFF or NetCDF into a Jupyter notebook in order to conduct further analysis or combine it with other Digital Earth Australia (DEA) data. In this example, we will demonstrate how to load in one or multiple GeoTIFF or NetCDF files originally exported to files from a Landsat-8 time-series into an xarray.Dataset for further analysis.

For advice on exporting raster data, refer to the Exporting GeoTIFFs notebook.

For more information on the xarray functions used:

xarray.open_rasterio (http://xarray.pydata.org/en/stable/generated/xarray.open_rasterio.html)
xarray.open_dataset (http://xarray.pydata.org/en/stable/generated/xarray.open_dataset.html)
xarray.open_mfdataset (http://xarray.pydata.org/en/stable/generated/xarray.open_mfdataset.html)

Description¶

This notebook shows how to open raster data from file using xarray’s built-in fuctions for handling GeoTIFF and NetCDF files:

Opening single raster files
- Opening a single GeoTIFF file
- Opening a single NetCDF file
Opening multiple raster files as an xarray.Dataset with a time dimension
- Opening multiple GeoTIFF files
- Opening multiple NetCDF files

Getting started¶

To run this example, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages¶

[1]:

%matplotlib inline

# import datacube
import sys
import glob
import xarray as xr

sys.path.append('../Scripts')
from dea_plotting import rgb
from dea_datahandling import paths_to_datetimeindex

/g/data/v10/public/modules/dea/20200526/lib/python3.6/site-packages/datacube/storage/masking.py:4: DeprecationWarning: datacube.storage.masking has moved to datacube.utils.masking
  category=DeprecationWarning)

Loading multiple files into a single xarray.Dataset¶

Geospatial time series data is commonly stored as multiple individual files with one time-step per file. These are difficult to use individually, so it can be useful to load multiple files into a single xarray.Dataset prior to analysis. This also allows us to analyse data in a format that is directly compatible with data directly loaded from the Datacube.

Multiple GeoTIFFs¶

To load multiple GeoTIFF files into a single xarray.Dataset, we first need to obtain a list of the files using the glob package. In the example below, we return a list of any files that match the pattern geotiff_*.tif:

[7]:

geotiff_list = glob.glob('../Supplementary_data/Opening_GeoTIFFs_NetCDFs/geotiff_*.tif')
geotiff_list

[7]:

['../Supplementary_data/Opening_GeoTIFFs_NetCDFs/geotiff_red_2018-01-19.tif',
 '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/geotiff_red_2018-01-03.tif',
 '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/geotiff_red_2018-03-24.tif']

We now read these files using xarray. To ensure each raster is labelled correctly with its time, we can use the helper function paths_to_datetimeindex() from dea_datahandling to extract time information from the file paths we obtained above. We then load and concatenate each dataset along the time dimension using xarray.open_rasterio(), convert the resulting xarray.DataArray to a xarray.Dataset, and give the variable a more useful name (red):

[8]:

# Create variable used for time axis
time_var = xr.Variable('time', paths_to_datetimeindex(geotiff_list,
                                                      string_slice=(12, -4)))

# Load in and concatenate all individual GeoTIFFs
geotiffs_da = xr.concat([xr.open_rasterio(i) for i in geotiff_list],
                        dim=time_var)

# Covert our xarray.DataArray into a xarray.Dataset
geotiffs_ds = geotiffs_da.to_dataset('band')

# Rename the variable to a more useful name
geotiffs_ds = geotiffs_ds.rename({1: 'red'})

# Print the output
print(geotiffs_ds)

<xarray.Dataset>
Dimensions:  (time: 3, x: 191, y: 212)
Coordinates:
  * y        (y) float64 -3.307e+06 -3.307e+06 ... -3.313e+06 -3.313e+06
  * x        (x) float64 2.053e+06 2.053e+06 2.053e+06 ... 2.058e+06 2.058e+06
  * time     (time) datetime64[ns] 2018-01-19 2018-01-03 2018-03-24
Data variables:
    red      (time, y, x) int16 902 1307 931 416 544 873 ... 244 285 307 251 247
Attributes:
    transform:      (30.0, 0.0, 2052720.0, 0.0, -30.0, -3306750.0)
    crs:            +init=epsg:3577
    res:            (30.0, 30.0)
    is_tiled:       1
    nodatavals:     (0.0,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

To verify the data was loaded correctly, we can plot it using xarray:

[9]:

geotiffs_ds.red.plot(col='time')

[9]:

<xarray.plot.facetgrid.FacetGrid at 0x7fb5b0084400>

../../_images/notebooks_Frequently_used_code_Opening_GeoTIFFs_NetCDFs_24_1.png

Multiple NetCDFs¶

The xarray.open_mfdataset() function can be used to easily load multiple NetCDFs into a single xarray.Dataset. First, we obtain file paths for the files we want to load using glob:

[10]:

netcdf_list = glob.glob('../Supplementary_data/Opening_GeoTIFFs_NetCDFs/netcdf_*.nc')
netcdf_list

[10]:

['../Supplementary_data/Opening_GeoTIFFs_NetCDFs/netcdf_2018-01-19.nc',
 '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/netcdf_2018-03-24.nc',
 '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/netcdf_2018-01-03.nc']

We then load in the NetCDF files using xarray.open_mfdataset. Because the NetCDF file format already contains time information for each dataset, we do not need to set up a time variable like in the previous GeoTIFF example.

[11]:

netcdfs_ds = xr.open_mfdataset(paths=netcdf_list, combine='by_coords')
netcdfs_ds

[11]:

xarray.Dataset

Dimensions:
- time: 3
- x: 191
- y: 212

Coordinates: (3)

y
(y)
float64
-3.307e+06 ... -3.313e+06
units :
metre
standard_name :
projection_y_coordinate
long_name :
y coordinate of projection
```
array([-3306765., -3306795., -3306825., ..., -3313035., -3313065., -3313095.])
```

x

(x)

float64

2.053e+06 2.053e+06 ... 2.058e+06

units :: metre
standard_name :: projection_x_coordinate
long_name :: x coordinate of projection

array([2052735., 2052765., 2052795., 2052825., 2052855., 2052885., 2052915.,
       2052945., 2052975., 2053005., 2053035., 2053065., 2053095., 2053125.,
       2053155., 2053185., 2053215., 2053245., 2053275., 2053305., 2053335.,
       2053365., 2053395., 2053425., 2053455., 2053485., 2053515., 2053545.,
       2053575., 2053605., 2053635., 2053665., 2053695., 2053725., 2053755.,
       2053785., 2053815., 2053845., 2053875., 2053905., 2053935., 2053965.,
       2053995., 2054025., 2054055., 2054085., 2054115., 2054145., 2054175.,
       2054205., 2054235., 2054265., 2054295., 2054325., 2054355., 2054385.,
       2054415., 2054445., 2054475., 2054505., 2054535., 2054565., 2054595.,
       2054625., 2054655., 2054685., 2054715., 2054745., 2054775., 2054805.,
       2054835., 2054865., 2054895., 2054925., 2054955., 2054985., 2055015.,
       2055045., 2055075., 2055105., 2055135., 2055165., 2055195., 2055225.,
       2055255., 2055285., 2055315., 2055345., 2055375., 2055405., 2055435.,
       2055465., 2055495., 2055525., 2055555., 2055585., 2055615., 2055645.,
       2055675., 2055705., 2055735., 2055765., 2055795., 2055825., 2055855.,
       2055885., 2055915., 2055945., 2055975., 2056005., 2056035., 2056065.,
       2056095., 2056125., 2056155., 2056185., 2056215., 2056245., 2056275.,
       2056305., 2056335., 2056365., 2056395., 2056425., 2056455., 2056485.,
       2056515., 2056545., 2056575., 2056605., 2056635., 2056665., 2056695.,
       2056725., 2056755., 2056785., 2056815., 2056845., 2056875., 2056905.,
       2056935., 2056965., 2056995., 2057025., 2057055., 2057085., 2057115.,
       2057145., 2057175., 2057205., 2057235., 2057265., 2057295., 2057325.,
       2057355., 2057385., 2057415., 2057445., 2057475., 2057505., 2057535.,
       2057565., 2057595., 2057625., 2057655., 2057685., 2057715., 2057745.,
       2057775., 2057805., 2057835., 2057865., 2057895., 2057925., 2057955.,
       2057985., 2058015., 2058045., 2058075., 2058105., 2058135., 2058165.,
       2058195., 2058225., 2058255., 2058285., 2058315., 2058345., 2058375.,
       2058405., 2058435.])

time

(time)

datetime64[ns]

2018-01-03T23:42:39 ... 2018-03-24T23:42:01

standard_name :: time
long_name :: Time, unix time-stamp
axis :: T

array(['2018-01-03T23:42:39.000000000', '2018-01-19T23:42:31.000000000',
       '2018-03-24T23:42:01.000000000'], dtype='datetime64[ns]')

Data variables: (4)

crs
(time)
int32
-2147483647 -2147483647 -2147483647
grid_mapping_name :
albers_conical_equal_area
standard_parallel :
[-18. -36.]
longitude_of_central_meridian :
132.0
latitude_of_projection_origin :
0.0
false_easting :
0.0
false_northing :
0.0
long_name :
GDA94 / Australian Albers
semi_major_axis :
6378137.0
semi_minor_axis :
6356752.314140356
inverse_flattening :
298.257222101
crs_wkt :
PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]
spatial_ref :
PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]
GeoTransform :
[ 2.05272e+06 3.00000e+01 0.00000e+00 -3.30675e+06 0.00000e+00 -3.00000e+01]
```
array([-2147483647, -2147483647, -2147483647], dtype=int32)
```

nbart_red

(time, y, x)

float32

dask.array<chunksize=(1, 212, 191), meta=np.ndarray>

grid_mapping :: crs
units :: 1





  
      Array  Chunk 
  
  
     Bytes  485.90 kB   161.97 kB 
     Shape  (3, 212, 191)   (1, 212, 191) 
     Count  12 Tasks  3 Chunks 
     Type  float32  numpy.ndarray

nbart_green

(time, y, x)

float32

dask.array<chunksize=(1, 212, 191), meta=np.ndarray>

grid_mapping :: crs
units :: 1





  
      Array  Chunk 
  
  
     Bytes  485.90 kB   161.97 kB 
     Shape  (3, 212, 191)   (1, 212, 191) 
     Count  12 Tasks  3 Chunks 
     Type  float32  numpy.ndarray

nbart_blue

(time, y, x)

float32

dask.array<chunksize=(1, 212, 191), meta=np.ndarray>

grid_mapping :: crs
units :: 1





  
      Array  Chunk 
  
  
     Bytes  485.90 kB   161.97 kB 
     Shape  (3, 212, 191)   (1, 212, 191) 
     Count  12 Tasks  3 Chunks 
     Type  float32  numpy.ndarray

Attributes: (11)
date_created :
2019-12-04T16:19:10.996947
Conventions :
CF-1.6, ACDD-1.3
history :
NetCDF-CF file created by datacube version '1.7+142.g7f8581cf' at 20191204.
geospatial_bounds :
POLYGON ((153.390117090279 -28.8513061571954,153.401115563742 -28.9072132711656,153.459739141404 -28.8986842614383,153.448711679615 -28.8427816411366,153.390117090279 -28.8513061571954))
geospatial_bounds_crs :
EPSG:4326
geospatial_lat_min :
-28.907213271165624
geospatial_lat_max :
-28.842781641136636
geospatial_lat_units :
degrees_north
geospatial_lon_min :
153.39011709027872
geospatial_lon_max :
153.45973914140396
geospatial_lon_units :
degrees_east

Because the NetCDF files we loaded using xarray contain multiple satellite bands (e.g. nbart_red, nbart_green, nbart_blue), we can plot the result in true colour:

[12]:

rgb(netcdfs_ds, col='time', percentile_stretch=(0.0, 0.9))

../../_images/notebooks_Frequently_used_code_Opening_GeoTIFFs_NetCDFs_30_0.png

Additional information¶

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on Github.

Last modified: June 2020

Tags¶

Browse all available tags on the DEA User Guide’s Tags Index

Tags: sandbox compatible, NCI compatible, dea_plotting, dea_datahandling, rgb, paths_to_datetimeindex, GeoTIFF, NetCDF, external data, xarray.open_rasterio, xarray.open_dataset, xarray.open_mfdataset

Opening GeoTIFF and NetCDF files with xarray ¶

Background¶

Description¶

Getting started¶

Load packages¶

Opening a single raster file¶

Define file paths¶

Opening a single GeoTIFF¶

Opening a single NetCDF¶

Loading multiple files into a single xarray.Dataset¶

Multiple GeoTIFFs¶

Multiple NetCDFs¶

Additional information¶

Tags¶