join_two_rasters_netcdf.py

import subprocess

import xarray as xr


# convert input tifs to netcdf
subprocess.check_call(['gdal_translate', '-of', 'NetCDF', 'biomass.tif', 'biomass.nc'])
subprocess.check_call(['gdal_translate', '-of', 'NetCDF', 'extent.tif', 'extent.nc'])

# rename bands to match data
subprocess.check_call(['ncrename', '-v', 'Band1,biomass', 'biomass.nc'])
subprocess.check_call(['ncrename', '-v', 'Band1,extent', 'extent.nc'])

# run a netcdf append command to add biomass data to the extent netcdf
subprocess.check_call(['ncks', '-A', 'biomass.nc', 'extent.nc'])

# open extent and then write it directly to CSV
# unfortunately there doesn't seem to be a fast CLI tool for this
ds = xr.open_dataset('extent.nc')
df = ds.to_dataframe().reset_index()
print df.head()

df.to_csv('netcdf_output.csv', index=False)