Code is publicly available at

Stitching to Save the Ice: Final Report

Greg Dray, Luke Alvoeiro and Tricia Nelsen


        The aim of our project is to stitch together images of the Arctic that are captured by the MODIS satellite for any given day. Additionally, in this combined image, we plan to flag pixels that most likely represent cloud cover. Images of the Arctic are essential to understanding the effect of the changing climate on the environment; the reflectivity of the ice forms a positive feedback system raising global temperatures. We chose this project because Tricia did related work this past summer and we hoped to assist her co-workers in their research by allowing them to perform more comprehensive analyses of sea ice floes.


        It is no secret that the Arctic is melting. The state of the Arctic Ocean is a strong indicator of global climate, and if the melt is perpetuated will lead to a warmer climate. The more open ocean present, the more incoming solar radiation is absorbed into the ocean, which directly impacts the rest of the global climate system (Perovich, et al., 2017). Each year the Arctic Ocean goes through stages of growth in the winter months and melt in the summer months, with the maximum extent usually two to three times larger than the minimum extent (Perovich, et al., 2017). Each year, the sea ice extent reaches its minimum around the end of September, and the maximum around March (Perovich, et al., 2017). Data has been recorded since 1979 on the extents of the Arctic sea ice. Since then, a clear decline has been seen in both the maximum and minimum extents of the ice. The lowest twelve extents have all occurred in the past twelve years, with this year’s extent tieing for the sixth lowest ice extent over the 40-year record (Gautier, 2018). The photo below shows the minimum extent for 2018 – 4.59 million square kilometers on September 19 – compared to the 1981-2010 September median extent, which is a metric used to compare annual sea ice extents (Viñas & Carlowicz, 2018).

Arctic Sea Ice Reaches 2018 Minimum

Similarly, the lowest four maximum ice extents have all occurred over the past four years, with 2015-2017 setting new record lows each year (Scott, 2018). These declining trends in Arctic sea ice are concerning and are crucial to understand.

In order to mitigate the increasing melt rate of the Arctic, the cause of the melt needs to be further researched. This is a challenging problem to solve due to the inaccessibility of sea ice. Some field work is being done to fully understand sea ice dynamics, however satellite imagery is a key component in being able to extract current conditions on a given day.

One tool that is being developed by Ph.D. candidate Nicholas Wright at Dartmouth College Thayer School of Engineering attempts to classify melt ponds on sea ice floes to assess the state of the ice and calculate the albedo (reflectivity) of a given floe. His algorithm is currently producing very accurate results, however, is only being applied to high-resolution images encompassing low spatial ranges. This means that, while the results are promising, large-scale conclusions about the state of the entire Arctic region cannot be made. There is a need for a product that can create an image of the entire Arctic on any given day to subsequently analyze the state of the sea ice.



        To create this product we decided to use MODIS (Moderate Resolution Imaging Spectroradiometer) data. MODIS is an instrument aboard two of NASA’s satellites that captures daily imagery of Earth’s surface (MODIS). Specifically, we are using product MOD09GA, a seven-band surface reflectance product available on a daily basis (Vermote, Kotchenova, & Ray, 2011). To extract information about sea ice, we are mainly concerned with reflectance within the visible light range: MODIS bands 1(red), 3 (blue), and 4 (green), at a 500-meter resolution.

In order to stitch the images together, we first needed to download the satellite images corresponding to the tiles and visible bands we desired. We went through the following libraries to perform this:

  1. MODIStsp: a package for R, was initially chosen as a source to automate the download of NASA MODIS imagery. We wrote a script to compare the edge coordinates of two images. This showed us that the images were already aligned correctly for perfect stitching. We wanted to see if there was a more challenging alternative.
  2. pyModis: a library for Python, we hoped, would offer us a different set of non-aligned images. After modifying aspects of our original script, we discovered that pyModis provided us with similar, aligned images. Using pyModis, we then wrote an algorithm that joins all the images in the correct order based on filenames we created.

Cloud Detection – First Approach

We attempted to use the analysis discussed in “Automated Detection of Clouds in Satellite Imagery” to inform our cloud cover detection algorithm (2009). The figure below shows the reflectance values of different elements seen aerially on Earth’s surface across the wavelength spectrum. The blue line is snow, which has a very similar reflectance to ice, and snow often covers the surface of sea ice. The red and green line are different types of clouds, while the pink line is forest canopy and the black line is soil. We used these reflectance values for clouds and snow at different wavelengths to create thresholds to flag a pixel as a cloud or not for each band that we analyzed from MODIS.

The table below shows the MOD09 surface reflectance band, the wavelengths that band detects, and the threshold we created to flag a pixel as a cloud.

MOD09 Surface Reflectance Band


Reflectance Range to Flag Pixel as Cloud

1 (red)



2 (near infrared)



3 (blue)



4 (green)

.55 - .57


Unfortunately, this algorithm did not work well. It seemed quite overzealous in what it labeled as clouds and included many areas that were dark and clearly not clouds. We believe this could be caused by two factors. First, the reflectances of the two types of clouds included here is very similar to that of both snow and forest canopy, so it might have been mistaking the two. Second, this method only considered two types of clouds, and it is possible that additional important types exist with different reflectances.

Cloud Detection – Second Approach

        Our next approach at using our image processing skills to create a cloud mask was to create a normalized difference snow index (NDSI). NDSI is a method that takes advantage of the fact that both clouds and snow have high reflectance values in visible wavelengths, but snow has low reflectance values in the shortwave infrared (SWIR) wavelength (MODIS’s 6th band) and clouds have higher reflectance values in SWIR To create this index we took the values in band 6 and subtracted them from band 4, MODIS’s green band. Then we divided this by the sum of the band values. We thresholded the mask as values greater than .4 would be considered cloud. Unfortunately, this did not work either, it detected almost nothing as clouds so we moved on to another approach.

Cloud Detection – Final Approach

        After our first two custom approaches both failed, we decided to switch to what we deemed to be a more reliable source. After some research, we discovered that a band called State_1km: Reflectance Data State actually contained two bits representing cloud cover information. (See Table 16 in the MODIS User Guide.) These bits are based on NASA’s own cloud detection algorithm, which is likely much more accurate than something we could devise. They indicate that a pixel is 1) clear, 2) cloudy, 3) mixed, or 4) not set. So as to not be too aggressive when labeling cloudy pixels, we did not label pixels that were marked as ‘mixed’ as ‘cloudy’ in our final output.

        The process involved in getting this data is straightforward. We call pyModis to download the additional band containing this cloud information. On each pixel in this band, we perform bitwise operations to remove all pixels except the ones containing cloud info. We then create a mask with 1 if a pixel is cloudy, 0 else. One minor complication is that this band is only at 1000m x 1000m resolution, and our visible bands are at 500m x 500m. We circumvent this issue by using numpy.repeat(), which essentially stretches our mask to be double in size, duplicating every value. The final step is simply to apply this mask to our stitched image. We mark cloudy pixels in our stitched image as [255, 0, 0], since pure red does not occur naturally.


Following this cloud detection and image stitching we used GDAL, specifically the gdalwarp, to reproject the image to a polar stereographic projection (EPSG 3413). This projection shows the Arctic north of 60°N, as if you were looking down on the north pole from above.


Our script consists of the following functions. These functions can be called as a whole, producing specific outputs, or called individually, to allow personalized processing of Arctic images:

This is an example of the three different tiff files we had for a single tile that were             combined into a single stacked file using combineChannels().


Since our product was intended to be used by researchers, we have tailored it accordingly; the program is designed to be run from the command line with the following arguments:

Above is an example of a fully stitched image of the Arctic. This is one of the 4 outputs that the user receives after running our Python script.

Our final output consists of 5 files:

  1. Stitched image
  2. Stitched and reprojected image
  3. Stitched and masked image
  4. Stitched, masked and reprojected image
  5. Cloud mask (optional)

The above image corresponds to the cloud masked stitched file. Areas marked as red correspond to regions where the MODIS algorithm detects cloud cover.

The above image corresponds to the fully stitched and reprojected image.

The above image corresponds to the fully stitched, reprojected, and masked image.


Planned schedule:

October 16

October 31

November 30

December 8

Tricia & Greg: Stitching strategies researched.

Luke: Python script for download complete.

Stitching strategy implemented without overlap. Begin prototyping for larger scale/automation.

Stitching accomplished automatically, dealing with overlap.

Code cleaned up and turned into final product with basic command line interaction.

Our actual schedule changed quite a bit from what we originally planned. We procrastinated quite a bit and some steps took much longer than we intended. Given our issue of actually obtaining data that had overlapping pixels, we spent a lot of time going back to square one to find new data. By November 10th we finally had tiles downloaded through Python and were able to process them into GeoTiffs. A few days later we were able to combine the bands and show the individual tiles. By the 26th we had successfully implemented the stitching of the entire Arctic. Finally, by December 6th, we had developed our cloud masking algorithm. Work performed after this date aimed to provide a easier to use product through command line options and cleaner code.


While working on this project, we have run into the following roadblocks:

Future Work

If we have time left over after implementing the cloud detection aspect of our solution, we intend to work on blending side-by-side images. This may be necessary as the satellite does not take all pictures of the Arctic at the same time; there is up to a 22-hour delay between two given pictures. The uniform blending of the stitched images is important, particularly since subsequence algorithms rely on similar images being outputted on a daily basis.

Additionally, we would explore additional ways to represent cloudy pixels. Although our current approach is functional, it is more ‘hard-coded’ than we would like. We might want to allow more flexibility for how these pixels are marked, and maybe even include an option where the marking is not visible but is simply another channel in the image.


Gautier, A. (2018, September 27). Arctic Sea Ice Extent Arrives at its Minimum. Arctic Sea Ice News & Analysis. Retrieved November 3, 2018, from

Jedlovec, G. (2009). Automated Detection of Clouds in Satellite Imagery. In Advances in Geoscience and Remote Sensing. InTechOpen. doi:10.5772/8326

MODIS. (n.d.). Retrieved October 5, 2018, from

Perovich, D., Meier, W., Tschudi, M., Farrell, S., Hendricks, S., Gerland, S., Haas, C., Krumpen, T., Polashenski, C., Ricker, R., Webster, M., 2017: Sea Ice [in Arctic Report Card 2017],

Scott, M. (2018, March 26). Arctic Sea Ice Extent at 2018 Winter Maximum was Second Smallest on Record. NOAA Climate News & Features. Retrieved December 3, 2018, from

Vermote, E. F., Kotchenova, S. Y., & Ray, J. P. (2011, February). MODIS Surface Reflectance User’s Guide(Rep.). Retrieved October 2, 2018, from NASA website:

Viñas, M., & Carlowicz, M. (2018, September 19). Arctic Sea Ice Reaches 2018 Minimum. NASA Earth Observatory Image of the Day. Retrieved December 3, 2018, from