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
andDEA 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 dimensionOpening 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)
Opening a single raster file¶
Define file paths¶
In the code below we define the locations of the GeoTIFF and NetCDF files that we will open. These files were originally exported from the GA Landsat 8 ga_ls8c_ard_3
product. GeoTIFF files contain a single satellite band (nbart_red
), while NetCDF files contain three satellite bands (nbart_red
, nbart_green
, nbart_blue
).
[2]:
geotiff_path = '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/geotiff_red_2018-01-03.tif'
netcdf_path = '../Supplementary_data/Opening_GeoTIFFs_NetCDFs/netcdf_2018-01-03.nc'
Opening a single GeoTIFF¶
To open a geotiff we use xarray.open_rasterio()
function which is built around the rasterio
Python package. When dealing with extremely large rasters, this function can be used to load data as a Dask array by providing a chunks
parameter (e.g. chunks={'x': 1000, 'y': 1000}
). This can be useful to reduce memory usage by only loading the specific portion of the raster you are interested in.
[3]:
# Open into an xarray.DataArray
geotiff_da = xr.open_rasterio(geotiff_path)
# Covert our xarray.DataArray into a xarray.Dataset
geotiff_ds = geotiff_da.to_dataset('band')
# Rename the variable to a more useful name
geotiff_ds = geotiff_ds.rename({1: 'red'})
We can plot the data to verify it loaded correctly:
[4]:
geotiff_ds.red.plot()
[4]:
<matplotlib.collections.QuadMesh at 0x7fb5b045d7f0>

Opening a single NetCDF¶
To open a NetCDF file we can use xarray.open_dataset()
. Similarly to the GeoTIFF example above, this function can also be used to open large NetCDF files as Dask arrays by providing a chunks
parameter (e.g. chunks={'x': 1000, 'y': 1000}
).
[5]:
# Open into an xarray.DataArray
netcdf_ds = xr.open_dataset(netcdf_path)
# We can use 'squeeze' to remove the single time dimension
netcdf_ds = netcdf_ds.squeeze('time')
netcdf_ds
[5]:
- x: 191
- y: 212
- time()datetime64[ns]2018-01-03T23:42:39
- standard_name :
- time
- long_name :
- Time, unix time-stamp
- axis :
- T
array('2018-01-03T23:42:39.000000000', dtype='datetime64[ns]')
- 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)float642.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.])
- crs()int32...
- 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, dtype=int32)
- nbart_red(y, x)float32...
- grid_mapping :
- crs
- units :
- 1
array([[ 411., 422., 404., ..., 933., 670., 544.], [ 533., 539., 530., ..., 812., 666., 751.], [ 567., 609., 788., ..., 907., 840., 829.], ..., [ 269., 308., 412., ..., 319., 306., 314.], [ 278., 286., 428., ..., 321., 322., 326.], [ 304., 560., 1339., ..., 332., 322., 312.]], dtype=float32)
- nbart_green(y, x)float32...
- grid_mapping :
- crs
- units :
- 1
array([[ 625., 640., 595., ..., 940., 734., 716.], [ 794., 788., 763., ..., 935., 843., 799.], [ 827., 865., 1028., ..., 883., 864., 934.], ..., [ 457., 496., 600., ..., 446., 430., 441.], [ 487., 471., 599., ..., 473., 457., 462.], [ 473., 733., 1501., ..., 497., 481., 450.]], dtype=float32)
- nbart_blue(y, x)float32...
- grid_mapping :
- crs
- units :
- 1
array([[ 435., 430., 421., ..., 656., 558., 497.], [ 487., 528., 546., ..., 660., 592., 556.], [ 532., 604., 759., ..., 608., 639., 710.], ..., [ 326., 348., 421., ..., 344., 323., 326.], [ 348., 348., 515., ..., 352., 345., 341.], [ 371., 576., 1175., ..., 363., 357., 339.]], dtype=float32)
- date_created :
- 2019-12-04T16:19:10.855359
- 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 file we loaded using xarray
contain multiple satellite bands (e.g. nbart_red
, nbart_green
, nbart_blue
), we can plot the result in true colour:
[6]:
rgb(netcdf_ds)

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>

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]:
- time: 3
- x: 191
- y: 212
- 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)float642.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]')
- 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)float32dask.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)float32dask.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)float32dask.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
- 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))

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