Loading Modis NDVI time series into PostGIS raster

One of the motivations behind the development of PostGIS raster is the need to effectively handle time series of large raster layers. This has always represented something of a challenge for desktop GIS, particularly once the layers become too large to fit into memory. Modis NDVI data is a good example.

NDVI2003

The data can be downloaded from the NASA reverb site

http://reverb.echo.nasa.gov/reverb

Using the reverb “shopping cart” provides a list of files that can be downloaded using wget with the -i option. I.e

wget -i data_url_script_2012-12-08_112521.txt

Nasa provides time series of NDVI at a range of temporal and spatial resolutions For work over a large region the 1km resolution is good enough.

The easiest way to handle the files is to use Nasa’s own MRT tool (Modis Resampling Tool), downloaded from the same site. This is Java based. It requires a system restart on Linux to register the path for using commands after installing.

Once all the files have been downloaded I filtered out the HDFs using the following shell command.

 mkdir hdfs
 ls|grep "hdf$"|xargs -I '{}'  mv {} hdfs

After entering the hdfs directory I used R to direct the mosaicing and resampling. I wanted to extract a list of unique dates from the files first. This involved splitting the file names and extracting parts of the string.

  setwd("/media/extdrive/Modis/hdfs")
  f<-function(x)unlist(strsplit(x,"\\."))[2]
  fnames<-dir(pattern="*hdf*")
  dts<-unlist(lapply(fnames,f))
  yrs<-substr(dts,2,5)
  datestring<-strptime(paste(yrs, substr(dts,6,8)), format="%Y %j")
  # f2<-function(x)substr(unlist(strsplit(x,"\\."))[5],1,7)
  # dts<-unlist(lapply(fnames,f2))
  lst<-unique(dts)
  date_lst<-unique(datestring)

Now, that done, carefully pasting into the MRT commands within a loop results in reprojected Geotiffs for each month. It does take a little effort to get right as MRT is quite sensitive to the format of it’s prm project files so cat has to be used with just the right spacing.

finaldest<-"/home/duncan/rasters/modis2/"
  for (i in 1:length(lst)){
  fls<-paste(fnames[dts%in%lst[i]],"\n",sep="")
  outfile<-"temp.hdf"
  cat(fls,file="temp.txt")
  command<-paste("mrtmosaic -i",getwd(),"/temp.txt -s 1 -o ",outfile,sep="")
  system(command)
  outfile<-paste(finaldest,date_lst[i],".tif",sep="")
  cat("    
  INPUT_FILENAME = temp.hdf
        
  OUTPUT_FILENAME = ",outfile," 
  
  RESAMPLING_TYPE = NEAREST_NEIGHBOR
     
  OUTPUT_PROJECTION_TYPE = GEO
  
  DATUM = WGS84 
  
  ",file="temp.prm")
  
  system("resample -p temp.prm")
  
  }

Now comes the more interesting bit.

Bborie very helpfully explained Postgresql partitioned tables to me. This is a great way to handle large amounts of identical data that you may want to query as a whole without ending up with a single unfeasibly large table. The Postgresql manual provides very clear guidance based around a similar example.
http://www.postgresql.org/docs/9.0/static/ddl-partitioning.html

So to use partitioning you first set up a parent table from which each subsequent table inherits its characteristics. This can be done directly from R if you want to.

library(RODBC)
con<-odbcConnect("modelling")
query<-"create schema modis2;"
odbcQuery(con,query)
query<-"CREATE TABLE modis2.ndvi (
    id serial NOT NULL PRIMARY KEY,
    rid         int not null,
    fdate         date not null,
    rast        raster
);"
odbcQuery(con,query)

Now each sub table needs to be set up as inheriting from the parent table with constraints. I used a table for each year whose values were set from a temporary table used during import. Again careful pasting in R is needed in order to add the appropriate elements to the commands and queries within a loop. The dates are first loaded to the temp file using the -F flag, which adds the file names. The previous code had given the files correctly formatted dates, so it was then just a case of stripping out the first characters and coercing them to dates.

The index is built on the sub-tables.

setwd(finaldest)
fls<-dir()
yrs<-substring(fls,1,4)
for(i in unique(yrs)){
 query<-paste("create table modis2.ndvi_",i,"(CHECK (fdate >= DATE '",i,"-01-01' AND fdate <= DATE '",i,"-12-31' )) inherits (modis2.ndvi);",sep="")
  print(query)
 odbcQuery(con,query)
  subfls<-fls[grep(i,fls)]
  for (f in subfls){
    command<-paste("raster2pgsql -s 4326 -d  -M -R /home/duncan/rasters/modis2/",f," -F -t 500x500 temp|psql -d models",sep="")
    print(command)
    system(command)
    query<-"alter table temp add column fdate date; update temp set fdate = substring(filename,1, 10)::date;"
    print(query)
    odbcQuery(con,query)
    query<-paste("insert into modis2.ndvi_",i," (rid,fdate,rast) select rid,fdate,rast from temp;",sep="")
    print(query)
    odbcQuery(con,query)
  }
 query<-paste("CREATE INDEX modis_ndvi-",i,"_spindex ON modis2.ndvi_",i," using gist(st_convexhull(rast));",sep="")
 print(query)
 odbcQuery(con,query)
}


In the next post I will show how these layers can be queried in PostGIS. Although it is still not possible to visualise PostGIS raster layers clearly in QGIS, the QGIS value tool can be a useful way of looking quickly at the original Geotiffs.

NDVI2

3 thoughts on “Loading Modis NDVI time series into PostGIS raster

    • Presumably you mean when that you see nothing when you try to visualise the results in QGIS. This is because the gdal driver still doesn’t seem to work with out of database rasters (last time I checked). However, the real point of having the rasters in the database is to allow you to work with the data they contain, rather than visualise them. I use Geoserver for that. You can visualise the rasters one way and extract information in another. It would be nice if this wasn’t necessary, but it is not a big problem.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s