NASA MODIS data

So something I recently worked on was looking at NASA satellite images. As is with the scientific method, I was originally curious about snow cover on land in a particular area, and wondered what kind of information could be learned from looking at satellite images. A quick search lead me to the MODIS data that NASA collects, and publishes.

With snow cover information, it allows you to form a view onto how much water inflow there might be in the future. But how does one quantify snow cover? From the point of view of a satellite, looking down on an area and capturing an image, you will get:

Now how do you determine what is snow cover?

Normalized Difference Snow Index (NDSI)

Luckily we have this numerical indicator that a smart person came up with, that gives us a ratio of spectral bands over an area, which in turn gives us a ratio of how much radiation is reflected in each band. Snow absorbs short wave infrared, and cloud does not, which is the key way of determining the difference between snow and cloud. This indicator is mainly used in looking at glacier movement, being able to accurate track the effects of global warming. NASA has a publication on it of course.

Again, NASA has a great tool that allows you to select your area of interest, and then it will find corresponding tiles for you. And this is the beginning of exploring GIS. I am not knowledgeable in GIS by any means, and there is definitely a lot of room for improvement.

The data

The format of the files are HDF, which is basically just a binary file format designed to organise large amounts of data. For the MODIS files, the data structure is that of raster images, which is a matrix structure of pixels. To deal with the HDF file formats, I use the h4toh5 converter tool, to change it into a format that I can use with Python hdf library, which mainly deal with the h5 format. After all the data is in the proper format for Python, we are able to begin.

Looking at the area of interest

The first step is to only look at pixels that we care about in each tile. If you are lucky, then you might only have to look at a single tile that contains all your information. But sometimes the area you want to look at happens to need more than a single tile. The first step in grabbing the pixels of interest is to convert the MODIS projection into longitude and latitude.

I do this with the pyproj package:

    sinu = pyproj.Proj("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext")
    wgs84 = pyproj.Proj("+init=EPSG:4326")
    lon, lat = pyproj.transform(sinu, wgs84, xv, yv)

If you have a KML file, say from Google Maps or Google Earth, then you can convert it into a numerical polygon using shapely package. Then essentially what I’ve done is loop over all the points in a single satellite tile(s) lon/lat and check if it falls within the shape of interest, and output a 2D numpy array that indicates whether or not that pixel is one we care about. This gives us a masking layer to apply over the satellite images

Extracting pixels of interest

With the mask layer is effectively a filter, letting us extract pixels that we care about from each satellite image. The task is then to grab the NDSI values from the tile for the pixel that falls within our area of interest.

The rest of the steps involve doing smoothing and clustering, and then finally some averaging to give an output NDSI factor for a particular area of interest.

Packages used

    import h5py
    import pyproj
    import keytree
    import numpy as np
    import psycopg2
    from scipy import interpolate
    from shapely.geometry import box, shape
    from astropy.convolution import convolve, Gaussian2DKernel
    from sklearn.cluster import KMeans, AgglomerativeClustering, MiniBatchKMeans
    from sklearn.preprocessing import scale
    
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from mpl_toolkits.basemap import Basemap

Some of these will require local distributions, I did this all in Linux, so apt-get got me all the packages I needed.

Improvements

Is this the best way of solving the task I had? I’m not totally sure, I learned about PostGIS and spatial indexing recently, and it seems like a pretty useful tool, although maybe not needed in my application. The slowest part of my Python script is when I was extracting and processing pixels from each tile, so some improvements could be made there.

All in all, it was a fun exercise and it’s pretty neat being able to play with map objects.