Downloading and processing MODIS with fast and free tools

Here’s a simple way of processing MODIS: downloading, extracting the bands from HDF as GEOTIFF, merging tiles and reprojecting. Area extraction, NDVI calculation, cloud filtering and more can be done in GRASS for example. It’s Linux based, but can easily used in Windows. My way might not be the most elegant, but it’s fast and works fine.

The first step is to select the tiles, the time frame and product at http://reverb.echo.nasa.gov/reverb/. I take MOD09Q1 (250 m, 8 day, collection 6 reflectance) as example. Once orderd, we get a TXT file containing all file names. We download all files using wget. After download i mass rename all files to remove the creation date (using e.g. total commander) so the files are named like: MOD09Q1.A2000065.h21v05.006.

Then I extract Band 1, Band 2 and the Quality flag as GeoTiff. “aus” is the folder were the TIFFS will end up. Find out the name of the bands using gdalinfo.

wget -i modis.txt
aus="/media/data/MODIS/tiffs/"
for i in *.hdf; do gdal_translate -of GTiff 'HDF4_EOS:EOS_GRID:"'$i'":MOD_Grid_250m_Surface_Reflectance:sur_refl_b01' $aus`basename $i .hdf`_REF_band1.tif; done
for i in *.hdf; do gdal_translate -of GTiff 'HDF4_EOS:EOS_GRID:"'$i'":MOD_Grid_250m_Surface_Reflectance:sur_refl_b02' $aus`basename $i .hdf`_REF_band2.tif; done
for i in *.hdf; do gdal_translate -of GTiff 'HDF4_EOS:EOS_GRID:"'$i'":MOD_Grid_250m_Surface_Reflectance:sur_refl_qc_250m' $aus`basename $i .hdf`_QA.tif; done

I did not find a clean way for merging, so i suggest my “dirty” way, it needs to be repeated for each year and band, this is 5 minutes of copying but you can use it again and again and again for the next years then….here’s an example for band 1 year 2000. If you use a 16 day product every second line will give an error, but that you can ignore.

# 2000
gdal_merge.py -o 2000001_band1.tif *2000001*
gdal_merge.py -o 2000009_band1.tif *2000009*
gdal_merge.py -o 2000017_band1.tif *2000017*
gdal_merge.py -o 2000025_band1.tif *2000025*
gdal_merge.py -o 2000033_band1.tif *2000033*
gdal_merge.py -o 2000041_band1.tif *2000041*
gdal_merge.py -o 2000049_band1.tif *2000049*
gdal_merge.py -o 2000057_band1.tif *2000057*
gdal_merge.py -o 2000065_band1.tif *2000065*
gdal_merge.py -o 2000073_band1.tif *2000073*
gdal_merge.py -o 2000081_band1.tif *2000081*
gdal_merge.py -o 2000089_band1.tif *2000089*
gdal_merge.py -o 2000097_band1.tif *2000097*
gdal_merge.py -o 2000105_band1.tif *2000105*
gdal_merge.py -o 2000113_band1.tif *2000113*
gdal_merge.py -o 2000121_band1.tif *2000121*
gdal_merge.py -o 2000129_band1.tif *2000129*
gdal_merge.py -o 2000137_band1.tif *2000137*
gdal_merge.py -o 2000145_band1.tif *2000145*
gdal_merge.py -o 2000153_band1.tif *2000153*
gdal_merge.py -o 2000161_band1.tif *2000161*
gdal_merge.py -o 2000169_band1.tif *2000169*
gdal_merge.py -o 2000177_band1.tif *2000177*
gdal_merge.py -o 2000185_band1.tif *2000185*
gdal_merge.py -o 2000193_band1.tif *2000193*
gdal_merge.py -o 2000201_band1.tif *2000201*
gdal_merge.py -o 2000209_band1.tif *2000209*
gdal_merge.py -o 2000217_band1.tif *2000217*
gdal_merge.py -o 2000225_band1.tif *2000225*
gdal_merge.py -o 2000233_band1.tif *2000233*
gdal_merge.py -o 2000241_band1.tif *2000241*
gdal_merge.py -o 2000249_band1.tif *2000249*
gdal_merge.py -o 2000257_band1.tif *2000257*
gdal_merge.py -o 2000265_band1.tif *2000265*
gdal_merge.py -o 2000273_band1.tif *2000273*
gdal_merge.py -o 2000281_band1.tif *2000281*
gdal_merge.py -o 2000289_band1.tif *2000289*
gdal_merge.py -o 2000297_band1.tif *2000297*
gdal_merge.py -o 2000305_band1.tif *2000305*
gdal_merge.py -o 2000313_band1.tif *2000313*
gdal_merge.py -o 2000321_band1.tif *2000321*
gdal_merge.py -o 2000329_band1.tif *2000329*
gdal_merge.py -o 2000337_band1.tif *2000337*
gdal_merge.py -o 2000345_band1.tif *2000345*
gdal_merge.py -o 2000353_band1.tif *2000353*
gdal_merge.py -o 2000361_band1.tif *2000361*

# We can remove the other files now 
rm MOD*
# Then reproject everything to a geographic system, maybe into a new folder:
mkdir lat
aus="/media/data/MODIS/tiffs/lat/"

for i in *_band1.tif; do gdalwarp -t_srs EPSG:4326 $i $aus`basename $i .tif`.tif; done
for i in *_band2.tif; do gdalwarp -t_srs EPSG:4326 $i $aus`basename $i .tif`.tif; done
for i in ???????_QA.tif; do gdalwarp -t_srs EPSG:4326 $i $aus`basename $i .tif`.tif; done

What we have now can already be used, but you can continue in GRASS to calculate NDVI and filter the data.

 

Recent woody vegetation trends in Sahel

Our new paper looks at recent dynamics in woody vegetation in Sahel and finds some interesting patterns which are mainly controlled by human population density.

Martin Brandt, Pierre Hiernaux, Kjeld Rasmussen, Cheikh Mbow, Laurent Kergoat, Torbern Tagesson, Yahaya Ibrahim, Abdoulaye Wele, Compton J. Tucker, Rasmus Fensholt. Assessing woody vegetation trends in Sahelian drylands using MODIS based seasonal metrics. Remote Sensing of Environment, 2016, 183, 215-225.

  • Woody cover trends are estimated for Sahel based on MODIS dry season metrics.
  • Interannual fluctuations in foliage density are attenuated to monitor woody plant trends.
  • Increases (decreases) are seen in areas of low (high) human population.
  • Recent decreases only partially offset a general post-drought increase in Sahelian woody cover.

Woody plants play a major role for the resilience of drylands and in peoples’ livelihoods. However, due to their scattered distribution, quantifying and monitoring woody cover over space and time is challenging. We develop a phenology driven model and train/validate MODIS (MCD43A4, 500 m) derived metrics with 178 ground observations from Niger, Senegal and Mali to estimate woody cover trends from 2000 to 2014 over the entire Sahel at 500 m scale.

Over the 15 year period we observed an average increase of 1.7 (± 5.0) woody cover (%) with large spatial differences: No clear change can be observed in densely populated areas (0.2 ± 4.2), whereas a positive change is seen in sparsely populated areas (2.1 ± 5.2). Woody cover is generally stable in cropland areas (0.9 ± 4.6), reflecting the protective management of parkland trees by the farmers. Positive changes are observed in savannas (2.5 ± 5.4) and woodland areas (3.9 ± 7.3).

The major pattern of woody cover change reveals strong increases in the sparsely populated Sahel zones of eastern Senegal, western Mali and central Chad, but a decreasing trend is observed in the densely populated western parts of Senegal, northern Nigeria, Sudan and southwestern Niger. This decrease is often local and limited to woodlands, being an indication of ongoing expansion of cultivated areas and selective logging.

We show that an overall positive trend is found in areas of low anthropogenic pressure demonstrating the potential of these ecosystems to provide services such as carbon storage, if not over-utilized. Taken together, our results provide an unprecedented synthesis of woody cover dynamics in the Sahel, and point to land use and human population density as important drivers, however only partially and locally offsetting a general post-drought increase.

graph_abstract

Assessing Future Vegetation Trends and Restoration Prospects in the Karst Regions of Southwest China

Our latest article is not located in the Sahel, however, the method to assess the future persistence of vegetation trends is highly interesting in the context of ecosystem stability and resistance. The article is open access and freely available.

Tong, Xiaowei; Wang, Kelin; Brandt, Martin; Yue, Yuemin; Liao, Chujie; Fensholt, Rasmus. 2016. “Assessing Future Vegetation Trends and Restoration Prospects in the Karst Regions of Southwest China.” Remote Sens. 8, no. 5: 357.

 

To alleviate the severe rocky desertification and improve the ecological conditions in Southwest China, the national and local Chinese governments have implemented a series of Ecological Restoration Projects since the late 1990s. In this context, remote sensing can be a valuable tool for conservation management by monitoring vegetation dynamics, projecting the persistence of vegetation trends and identifying areas of interest for upcoming restoration measures.

In this study, we use MODIS satellite time series (2001–2013) and the Hurst exponent to classify the study area (Guizhou and Guangxi Provinces) according to the persistence of future vegetation trends (positive, anti-persistent positive, negative, anti-persistent negative, stable or uncertain). The persistence of trends is interrelated with terrain conditions (elevation and slope angle) and results in an index providing information on the restoration prospects and associated uncertainty of different terrain classes found in the study area.

The results show that 69% of the observed trends are persistent beyond 2013, with 57% being stable, 10% positive, 5% anti-persistent positive, 3% negative, 1% anti-persistent negative and 24% uncertain. Most negative development is found in areas of high anthropogenic influence (low elevation and slope), as compared to areas of rough terrain. We further show that the uncertainty increases with the elevation and slope angle, and areas characterized by both high elevation and slope angle need special attention to prevent degradation. Whereas areas with a low elevation and slope angle appear to be less susceptible and relevant for restoration efforts (also having a high uncertainty), we identify large areas of medium elevation and slope where positive future trends are likely to happen if adequate measures are utilized.

The proposed framework of this analysis has been proven to work well for assessing restoration prospects in the study area, and due to the generic design, the method is expected to be applicable for other areas of complex landscapes in the world to explore future trends of vegetation.

graphical_abstract

The Hurst exponent, a measure of the persistence of a time series, can be easily calculated on a raster time series in R:

library(raster)
library(rgdal)
library(pracma)

# set working directory with the raster files, here in TIF format, 
# and load the files in a rasterbrick

setwd("/media/2016_Xiaowei/")
gsn = brick(list.files(pattern='*.tif'))
# test on the average of the study area

g=cellStats(gsn, stat='mean')
gsn.ts = ts(g, start=c(2001,1), end=c(2013,1), frequency=1)
plot(gsn.ts)

h=hurstexp(gsn.ts)$Hs
h

# run the function on the rasterbrick gsn

fun=function(x) { 
  v=as.vector(x)
  if (is.na(v[1])){ NA } else
  gsn.ts = ts(v, start=c(2001,1), end=c(2013,1), frequency=1)
  x=hurstexp(gsn.ts, display=F)$Hs 
}
h <- calc(gsn, fun)
plot(h)

Renaming, converting, clipping: script based raster time series processing

When working with remotely sensed time series data (e.g. MODIS, AVHRR, GIMMS, etc.), bulk processing can save a lot of time. Using a terminal in a Linux environment, simple scripts can process thousends of files in short time. Here are some basic hints on how to start, gdal has to be installed. For Windows users or users who prefer an interface, QGIS provides a graphical interface for most of the gdal modules. Also renaming can be done more user friendly.

Renaming:

Original data-names can be very long annoying. So first we rename all files to numbers from 1 to XX (like 1, 2, 3…). To do so, navigate to the folder with the files, and these lines will rename all files originally starting with MOD (like MOD13Q1203837383922837284.hdf):

for fname in MOD* 
do
mv $fname `printf $x`
x=$(($x+1))
done

Then the following lines will create a script that renames all files numbered from 1 to 60. the final files will include year and date. So we then have data named MODIS_2001_01 until MODIS_2005_12. This can be adopted to every name and numbering. The script is then executed.

for ((i=1;$i<=60;i++)) do echo "mv $i MODIS_$(( ($i-1) / 12 + 2001 ))_$(( ($i-1) % 12 + 1 )).tif" ; done >rename.sh
sed -i 's/_1/_01/g' rename.sh
sed -i 's/_2/_02/g' rename.sh
sed -i 's/_3/_03/g' rename.sh
sed -i 's/_4/_04/g' rename.sh
sed -i 's/_5/_05/g' rename.sh
sed -i 's/_6/_06/g' rename.sh
sed -i 's/_7/_07/g' rename.sh
sed -i 's/_8/_08/g' rename.sh
sed -i 's/_9/_09/g' rename.sh
chmod 0755 rename.sh
./rename.sh

Convert:

This works fine to convert all kind of formats, also very useful to change HDF files (like MODIS, VGT…) to GeoTiff. This example converts GeoTiffs to ENVI format.

# define in and out folder and output format

in="/home/martin/SPOT_VGT/NDV/"
out="/home/martin/SPOT_VGT/NDV/converted/"
format="ENVI"

# convert all Tiffs to ENVI format, use float32 as output here

cd $in
for i in *.tif; do gdal_translate -of $format -ot Float32 $i $out`basename $i`; done

Clip:

# define in and out folder

in="/home/martin/SPOT_VGT/NDV/"
out="/home/martin/SPOT_VGT/NDV/clipped/"

# define the window which should be clipped

range="-17.6 16.25 -2.72 13.16"

# do the clipping for all Tiffs

cd $in 
for i in *.tif; do gdal_translate -projwin $range $i $out`basename $i`; done

For more information, see http://www.gdal.org/

Simple time series analysis with GIMMS NDVI and R

GIMMS NDVI time series for a selected pixel

GIMMS NDVI time series for a selected pixel

Time series analysis with satellite derived greenness indexes (e.g. NDVI) is a powerfull tool to assess environmental processes. AVHRR, MODIS and SPOT VGT provide global and daily imagery. Creating some plots is a simple task, and here is a rough start how it is done with GIMMS NDVI. All we need is the free software R.

I assume we already have the NDVI images in a folder, so the first 3 steps can be skipped, if not, they can be downloaded with a simple command:

wget -r ftp://ftp.glcf.umd.edu/glcf/GIMMS/Albers/Africa/

now unzip it:

find . -name "*.gz" -exec unp {} \;

As it comes in Albers projection, I reproject the Tiffs to a Lat Long WGS84 projection:

for in in *tif; do gdalwarp -s_srs EPSG:9001 -t_srs EPSG:4326 -r near -of GTiff $i /home/martin/Dissertation/GIMMS/GIMMS_lat/`basename $i`; done

I delete all 1981 images, as the year is not complete.

next, we start R, load the raster package and set the directory with the GIMMS Geotiffs:

library(raster)
library(rgdal)
setwd("~/GIMMS")

now load the Tiffs as a raster brick:

sg = stack(list.files(pattern='*.tif'))
gimmsb = brick(sg)
rm(sg)

now let’s chose a pixel and create a time series object. It can either be done by Row and Col:

i=20 # row
j=23 # col
gimms = gimmsb[i,j]
gimms.ts = ts(gimms, start=1982, end=c(2006,24), frequency=24)

or by coordinates

xy <- cbind(-3.6249,14.3844) 
sp <- SpatialPoints(xy) 
data <- extract(gimmsb, sp)
data=as.vector(data)
gimms.ts = ts(data, start=1982, end=c(2006,24), frequency=24)

now the time series can be used for further analysis or plotted:

plot(gimms.ts)

A regression line and a LOESS smoothing can be added:

gimms.time = time(gimms.ts) 
plot(gimms.ts)
abline(reg=lm(gimms.ts ~ gimms.time))
lines(lowess(gimms.ts, f=0.3, iter=5), col="blue")

some playing with boxplots and annual values:

layout(1:2)
plot(aggregate(gimms.ts)); 
boxplot(gimms.ts ~ cycle(gimms.ts))
GIMMS time series analysis

GIMMS time series analysis

averaged over the sudy area:

g=cellStats(gimmsb, stat='mean')
g.ts = ts(g, start=1982, end=c(2006,24), frequency=24)
plot(g.ts)

There are great possibilities for time series analysis in R, e.g. the series can be smoothed, decomposed…etc. The zoo package provides further functions and can handle irregular series. See also the lattice package for better visual presentation. How pixelwise analyis is conducted for time series of rasters, I’ll write in another post.

Detecting environmental change using time series, high resolution imagery and field work – a case study in the Sahel of Mali

Climatic changes and population pressure have caused major environmental change in the Sahel during the last fifty years. Many studies use coarse resolution NDVI time series such as GIMMS to detect environmental trends; however explanations for these trends remain largely unknown.

map

We suggest a five-step methodology for the validation of trends with a case study on the Dogon Plateau, Mali. The first step is to monitor long-term trends with coarse scale time series. Instead of GIMMS, we use a combination of LTDR (derived from AVHRR) and SPOT VGT NDVI data, covering the period from 1982-2010 with a temporal resolution of 10 days and a spatial resolution of 5 km.

Areas with significant trends are further analysed in a second step. Here we use a decomposed MODIS time series with a spatial resolution of 250 m to discover details of the blue spot i9n Figure 1. Due to the large scaled MODIS dataset, trends can be identified at a local scale / village level, see Figure 2.

Fig. 1: LTDR-SPOT showing spatial trends of NDVI. Spatial variations can be observed at a scale of 5.6 km (here the resolution is interpolated to 1 km). South of Fiko the large blue area stands out. This seems to be an area which does not show greening trends after the droughts in the beginning of the 80s.

Fig. 1: LTDR-SPOT showing spatial trends of NDVI. Spatial variations can be observed at a scale of 5.6 km (here the resolution is interpolated to 1 km). South of Fiko the large blue area stands out. This seems to be an area which does not show greening trends after the droughts in the beginning of the 80s.

Fig. 2: The map corresponds to the area of the rectangle in Fig. 1 and shows significant trends of MODIS time series since 2000. At a resolution of 250 m, the spatial patterns are far more diverse and variations within small areas can be detected. This demonstrates that the LTDR-SPOT trends merge many processes into single pixels. Thus further steps must be taken to explain local variations.

Fig. 2: The map corresponds to the area of the rectangle in Fig. 1 and shows significant trends of MODIS time series since 2000. At a resolution of 250 m, the spatial patterns are far more diverse and variations within small areas can be detected. This demonstrates that the LTDR-SPOT trends merge many processes into single pixels. Thus further steps must be taken to explain local variations.

Using very high resolution imagery (e.g. SPOT, Quickbird) areas of interest can be compared with pre-drought Corona-imagery from 1967. This offers a detailed overview of the environmental change at tree-level. A comparison of high resolution imagery with the Corona images show major land use changes over the past fifty years. What used to be dense bush cover has partially been converted to farmer managed agro-forestry and a significant proportion is now degraded land. Furthermore, an increase of tree cover on the fields can be detected. These different trends can also be observed in figures 3 and 4.

Fig. 3:  Left: Quickbird 2007 (Google Earth) Right: Corona 1967

Fig. 3: Left: Quickbird 2007 (Google Earth) Right: Corona 1967

Fig 4:  Left: Quickbird 2007 (Google Earth) Right: Corona 1967

Fig 4: Left: Quickbird 2007 (Google Earth) Right: Corona 1967

Yet many explanations for the changes identified remain unclear.

On-site field work provides information on the land use systems, vegetation composition and the current environmental condition. An initial field trip validated the suspected soil erosion and ongoing loss of trees and shrubs outside the fields used for farming purposes. On the fields surrounding the village many useful trees of all ages were identified. Still many explanations for change can only be speculated and hypothesized.

Fig. 5: There is a huge difference of farming areas (second row) and grazing areas (first row) which are adjacent (images taken in Nov. 2011).

Fig. 5: There is a huge difference of farming areas (second row) and grazing areas (first row) which are adjacent (images taken in Nov. 2011).

Fig. 6: On the one hand trees are protected on farmer's fields, on the other hand the „brousse / forêt“ is exploited for firewood (images taken in Nov. 2011).

Fig. 6: On the one hand trees are protected on farmer’s fields, on the other hand the „brousse / forêt“ is exploited for firewood (images taken in Nov. 2011).

Still many explanations for change can only be speculated and hypothesized. For this reason, interviews with the local population are vital for providing missing details.

Interviews with local people showed that good farmer-management using traditional methods, without outside-influence of projects, led to an increase of tree cover on the fields and healthy environmental conditions.

The land outside of the current farming area is highly degraded, which locals explain by the following points:

  • the extreme droughts in the 1970s and 1980s,
  • lack of rain in the past 30 years,
  • lack of protection by farmers,
  • legal and illegal felling by inhabitants of provincial towns in the region,
  • increased livestock numbers put pressure on the soil and vegetation.

Due to the declining vegetation cover and supported by the unfavourable morphology, the susceptibility to soil erosion increases. Many useful trees and shrubs have become rare or disappeared in these areas (e.g. Butyrospermum parkii, Crataeva adansonii, Combretum micranthum, Piliostigma reticulatum, Pterocarpus lucens, Sclerocarya birrea, etc).

Fig. 6: On the one hand trees are protected on farmer's fields, on the other hand the „brousse / forêt“ is exploited for firewood (images taken in Nov. 2011).

Fig. 6: On the one hand trees are protected on farmer’s fields, on the other hand the „brousse / forêt“ is exploited for firewood (images taken in Nov. 2011).

Fig. 7: Interviewing a farmer from Djamnati (image taken in Nov. 2011).

Fig. 7: Interviewing a farmer from Djamnati (image taken in Nov. 2011).

This example demonstrates the importance of land use and how an integrative and qualitative approach as well as input of local inhabitants expands knowledge and understanding of environmental change in the Sahel. Greening and degradation have many reasons which need to be varified by field work. Our example demonstrates, that climatic factors are important drivers of environmental changes. But land use concepts lead to oppositional results in vegetation development and therefore heterogenous landscape patterns.

see the poster:

egu_poster_small

see: Brandt, M., Samimi, C., Romankiewicz, C. & R. Spiekermann (2012): Detecting environmental change using time series, high resolution imagery and field work – a case study in the Sahel of Mali. Geophysical Research Abstracts, Vol. 14, EGU2012-10583, EGU General Assembly 2012.