In [1]:
%store -r
In [2]:
# Import libraries
import json
import os
import pathlib
import shutil
from glob import glob

import earthpy.api.appeears as eaapp
import earthpy
import geopandas as gpd
import hvplot.pandas
import hvplot.xarray
import pandas as pd
import rioxarray as rxr
import xarray as xr
In [3]:
# Get a sorted list of NDVI tif file paths
ndvi_paths = sorted(list(project.project_dir.rglob(
    'MOD13Q1.061__250m_16_days_NDVI*.tif')))

# Display the first and last three files paths to check the pattern
ndvi_paths[:3], ndvi_paths[-3:]
Out[3]:
([PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001145000000_aid0001.tif'),
  PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001161000000_aid0001.tif'),
  PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2001177000000_aid0001.tif')],
 [PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022209000000_aid0001.tif'),
  PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022225000000_aid0001.tif'),
  PosixPath('/workspaces/data/gc-tubarjal/tubarjal-ndvi/MOD13Q1.061_2001137_to_2022244/MOD13Q1.061__250m_16_days_NDVI_doy2022241000000_aid0001.tif')])
In [ ]:
doy_start = -25
doy_end = -18
# Loop through each NDVI image
ndvi_das = []
for ndvi_path in ndvi_paths:
    # Get date from file name
    # print(ndvi_path)
    date = pd.to_datetime(ndvi_path.name[doy_start:doy_end], format='%Y%j')
    # print(date)
    
    # Open dataset
    da = rxr.open_rasterio(ndvi_path,mask_and_scale=True).squeeze()
    # display(da)
    # da.plot()
    
    # Add date dimension and clean up metadata
    da = da.assign_coords({'date': date})
    da = da.expand_dims({'date': 1})
    da.name = 'NDVI'
    # display(da)
    # da.plot()
    
    # Prepare for concatenation
    ndvi_das.append(da)     
len(ndvi_das )

Combine Rasters¶

Next, stack your arrays by date into a time series using the xr.combine_by_coords() function. You will have to tell it which dimension you want to stack your data in.

In [ ]:
# Combine NDVI images from all dates
annual_ndvi_da = ( 
    xr.combine_by_coords(ndvi_das)
    .resample(dict(date = 'YS'))
    .mean()
)
annual_ndvi_da
Out[ ]:
<xarray.Dataset> Size: 1MB
Dimensions:      (date: 22, y: 85, x: 165)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 1kB 37.87 37.88 37.88 37.88 ... 38.21 38.21 38.22
  * y            (y) float64 680B 30.33 30.33 30.33 30.33 ... 30.16 30.16 30.16
    spatial_ref  int64 8B 0
  * date         (date) datetime64[ns] 176B 2001-01-01 2002-01-01 ... 2022-01-01
Data variables:
    NDVI         (date, y, x) float32 1MB 0.1336 0.1343 0.1339 ... 0.2285 0.1965
xarray.Dataset
    • date: 22
    • y: 85
    • x: 165
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      37.87 37.88 37.88 ... 38.21 38.22
      array([37.873958, 37.876042, 37.878125, 37.880208, 37.882292, 37.884375,
       37.886458, 37.888542, 37.890625, 37.892708, 37.894792, 37.896875,
       37.898958, 37.901042, 37.903125, 37.905208, 37.907292, 37.909375,
       37.911458, 37.913542, 37.915625, 37.917708, 37.919792, 37.921875,
       37.923958, 37.926042, 37.928125, 37.930208, 37.932292, 37.934375,
       37.936458, 37.938542, 37.940625, 37.942708, 37.944792, 37.946875,
       37.948958, 37.951042, 37.953125, 37.955208, 37.957292, 37.959375,
       37.961458, 37.963542, 37.965625, 37.967708, 37.969792, 37.971875,
       37.973958, 37.976042, 37.978125, 37.980208, 37.982292, 37.984375,
       37.986458, 37.988542, 37.990625, 37.992708, 37.994792, 37.996875,
       37.998958, 38.001042, 38.003125, 38.005208, 38.007292, 38.009375,
       38.011458, 38.013542, 38.015625, 38.017708, 38.019792, 38.021875,
       38.023958, 38.026042, 38.028125, 38.030208, 38.032292, 38.034375,
       38.036458, 38.038542, 38.040625, 38.042708, 38.044792, 38.046875,
       38.048958, 38.051042, 38.053125, 38.055208, 38.057292, 38.059375,
       38.061458, 38.063542, 38.065625, 38.067708, 38.069792, 38.071875,
       38.073958, 38.076042, 38.078125, 38.080208, 38.082292, 38.084375,
       38.086458, 38.088542, 38.090625, 38.092708, 38.094792, 38.096875,
       38.098958, 38.101042, 38.103125, 38.105208, 38.107292, 38.109375,
       38.111458, 38.113542, 38.115625, 38.117708, 38.119792, 38.121875,
       38.123958, 38.126042, 38.128125, 38.130208, 38.132292, 38.134375,
       38.136458, 38.138542, 38.140625, 38.142708, 38.144792, 38.146875,
       38.148958, 38.151042, 38.153125, 38.155208, 38.157292, 38.159375,
       38.161458, 38.163542, 38.165625, 38.167708, 38.169792, 38.171875,
       38.173958, 38.176042, 38.178125, 38.180208, 38.182292, 38.184375,
       38.186458, 38.188542, 38.190625, 38.192708, 38.194792, 38.196875,
       38.198958, 38.201042, 38.203125, 38.205208, 38.207292, 38.209375,
       38.211458, 38.213542, 38.215625])
    • y
      (y)
      float64
      30.33 30.33 30.33 ... 30.16 30.16
      array([30.332292, 30.330208, 30.328125, 30.326042, 30.323958, 30.321875,
       30.319792, 30.317708, 30.315625, 30.313542, 30.311458, 30.309375,
       30.307292, 30.305208, 30.303125, 30.301042, 30.298958, 30.296875,
       30.294792, 30.292708, 30.290625, 30.288542, 30.286458, 30.284375,
       30.282292, 30.280208, 30.278125, 30.276042, 30.273958, 30.271875,
       30.269792, 30.267708, 30.265625, 30.263542, 30.261458, 30.259375,
       30.257292, 30.255208, 30.253125, 30.251042, 30.248958, 30.246875,
       30.244792, 30.242708, 30.240625, 30.238542, 30.236458, 30.234375,
       30.232292, 30.230208, 30.228125, 30.226042, 30.223958, 30.221875,
       30.219792, 30.217708, 30.215625, 30.213542, 30.211458, 30.209375,
       30.207292, 30.205208, 30.203125, 30.201042, 30.198958, 30.196875,
       30.194792, 30.192708, 30.190625, 30.188542, 30.186458, 30.184375,
       30.182292, 30.180208, 30.178125, 30.176042, 30.173958, 30.171875,
       30.169792, 30.167708, 30.165625, 30.163542, 30.161458, 30.159375,
       30.157292])
    • spatial_ref
      ()
      int64
      0
      crs_wkt :
      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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      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
      grid_mapping_name :
      latitude_longitude
      spatial_ref :
      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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
      GeoTransform :
      37.87291666327381 0.0020833333331466974 0.0 30.333333330615915 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2001-01-01 ... 2022-01-01
      array(['2001-01-01T00:00:00.000000000', '2002-01-01T00:00:00.000000000',
       '2003-01-01T00:00:00.000000000', '2004-01-01T00:00:00.000000000',
       '2005-01-01T00:00:00.000000000', '2006-01-01T00:00:00.000000000',
       '2007-01-01T00:00:00.000000000', '2008-01-01T00:00:00.000000000',
       '2009-01-01T00:00:00.000000000', '2010-01-01T00:00:00.000000000',
       '2011-01-01T00:00:00.000000000', '2012-01-01T00:00:00.000000000',
       '2013-01-01T00:00:00.000000000', '2014-01-01T00:00:00.000000000',
       '2015-01-01T00:00:00.000000000', '2016-01-01T00:00:00.000000000',
       '2017-01-01T00:00:00.000000000', '2018-01-01T00:00:00.000000000',
       '2019-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000',
       '2021-01-01T00:00:00.000000000', '2022-01-01T00:00:00.000000000'],
      dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.1336 0.1343 ... 0.2285 0.1965
      units :
      NDVI
      AREA_OR_POINT :
      Area
      array([[[0.13361427, 0.1343    , 0.13387144, ..., 0.11654287,
         0.24504288, 0.3828    ],
        [0.13132857, 0.13165714, 0.13165714, ..., 0.21975715,
         0.21975715, 0.2943714 ],
        [0.12894285, 0.12892857, 0.12872855, ..., 0.18125713,
         0.19122855, 0.17778572],
        ...,
        [0.13184285, 0.13035713, 0.13035713, ..., 0.11554285,
         0.1148857 , 0.12895714],
        [0.13001429, 0.12939999, 0.1319    , ..., 0.11432857,
         0.12138572, 0.15157142],
        [0.1294    , 0.13069999, 0.13078572, ..., 0.10821428,
         0.12915714, 0.15702857]],

       [[0.13472857, 0.13417141, 0.1334    , ..., 0.12508571,
         0.32601425, 0.5295715 ],
        [0.13525714, 0.13511428, 0.13511428, ..., 0.2071    ,
         0.2071    , 0.3606143 ],
        [0.1356    , 0.13579999, 0.13705714, ..., 0.15152857,
         0.1651    , 0.18801428],
...
        [0.13257143, 0.12988572, 0.12988572, ..., 0.19052856,
         0.1374    , 0.21858573],
        [0.1270857 , 0.131     , 0.1334    , ..., 0.3250143 ,
         0.20534284, 0.14534284],
        [0.13012856, 0.13108571, 0.13509999, ..., 0.40951428,
         0.19429998, 0.1592    ]],

       [[0.13744286, 0.13808571, 0.1375    , ..., 0.1529857 ,
         0.12371429, 0.14465714],
        [0.13377143, 0.13592856, 0.13592856, ..., 0.14617144,
         0.14617144, 0.14565715],
        [0.13457142, 0.1369857 , 0.13571428, ..., 0.16374286,
         0.16337143, 0.13561428],
        ...,
        [0.13112856, 0.13147143, 0.13147143, ..., 0.2268857 ,
         0.14514287, 0.18932857],
        [0.12848571, 0.12875713, 0.13121428, ..., 0.47364283,
         0.26062855, 0.18037143],
        [0.12888572, 0.131     , 0.13152857, ..., 0.49725714,
         0.22854283, 0.19648571]]], shape=(22, 85, 165), dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([37.873958329940386,  37.87604166327353,  37.87812499660668,
       37.880208329939826,  37.88229166327297,  37.88437499660612,
       37.886458329939266,  37.88854166327241,  37.89062499660556,
       37.892708329938706,
       ...
       38.196874996578124,  38.19895832991127,  38.20104166324442,
       38.203124996577564,  38.20520832991071,  38.20729166324386,
       38.209374996577004,  38.21145832991015,   38.2135416632433,
       38.215624996576445],
      dtype='float64', name='x', length=165))
    • y
      PandasIndex
      PandasIndex(Index([ 30.33229166394934, 30.330208330616195, 30.328124997283048,
         30.3260416639499, 30.323958330616755, 30.321874997283608,
        30.31979166395046, 30.317708330617315, 30.315624997284168,
        30.31354166395102, 30.311458330617874, 30.309374997284728,
        30.30729166395158, 30.305208330618434, 30.303124997285288,
        30.30104166395214, 30.298958330618994, 30.296874997285848,
         30.2947916639527, 30.292708330619554, 30.290624997286407,
        30.28854166395326, 30.286458330620114, 30.284374997286967,
        30.28229166395382, 30.280208330620674, 30.278124997287527,
        30.27604166395438, 30.273958330621234, 30.271874997288087,
        30.26979166395494, 30.267708330621794, 30.265624997288647,
         30.2635416639555, 30.261458330622354, 30.259374997289207,
        30.25729166395606, 30.255208330622914, 30.253124997289767,
        30.25104166395662, 30.248958330623474, 30.246874997290327,
        30.24479166395718, 30.242708330624033, 30.240624997290887,
        30.23854166395774, 30.236458330624593, 30.234374997291447,
         30.2322916639583, 30.230208330625153, 30.228124997292007,
        30.22604166395886, 30.223958330625713, 30.221874997292566,
        30.21979166395942, 30.217708330626273, 30.215624997293126,
        30.21354166395998, 30.211458330626833, 30.209374997293686,
        30.20729166396054, 30.205208330627393, 30.203124997294246,
         30.2010416639611, 30.198958330627953, 30.196874997294806,
        30.19479166396166, 30.192708330628513, 30.190624997295366,
        30.18854166396222, 30.186458330629073, 30.184374997295926,
        30.18229166396278, 30.180208330629632, 30.178124997296486,
        30.17604166396334, 30.173958330630192, 30.171874997297046,
         30.1697916639639, 30.167708330630752, 30.165624997297606,
        30.16354166396446, 30.161458330631312, 30.159374997298166,
        30.15729166396502],
      dtype='float64', name='y'))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2001-01-01', '2002-01-01', '2003-01-01', '2004-01-01',
               '2005-01-01', '2006-01-01', '2007-01-01', '2008-01-01',
               '2009-01-01', '2010-01-01', '2011-01-01', '2012-01-01',
               '2013-01-01', '2014-01-01', '2015-01-01', '2016-01-01',
               '2017-01-01', '2018-01-01', '2019-01-01', '2020-01-01',
               '2021-01-01', '2022-01-01'],
              dtype='datetime64[ns]', name='date', freq='YS-JAN'))
In [ ]:
%store annual_ndvi_da

Stored 'annual_ndvi_da' (Dataset)

Finally, be sure to Restart and Run all to make sure your notebook works all the way through!

In [ ]:
%%capture
%%bash
#Coverting to HTML
jupyter nbconvert vegetation-02-wrangle.ipynb --to html