create capstone project

marwahaha commented 7 years ago

High-level goals:

Make Account: https://urs.earthdata.nasa.gov/users/new
pip install pymodis Problems? Try conda install gdal and try again Problems? Download the wheel here http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal
(Maybe skip this) Download, look at scripts (clone that repository.... or look at it?
Download EarthData
Talk about HDF formats
Mosaic, convert Hdf to GeoTiff
Talk about Open data formats, like GeoTiff
Open it in QGIS
Download shapefile, cut the shapefile

Initial commands that work...

python pyModis/scripts/modis_download.py -I -f 2012-12-05 -O -t h28v05,h29v05,h28v04 .
python pyModis/scripts/modis_mosaic.py -s "1" -o mosaic -v listfileMOD11A1.005.txt
python pyModis/scripts/modis_convert.py -v -s "( 1 )" -o final -e 4326 A2012340_mosaic_LST_Day_1km.vrt

http://www.wikience.org/tutorials/open-geotiff-in-quantum-gis/

marwahaha commented 7 years ago

And if you want to figure out shapefile cutting, you can use ogr2ogr -clipsrc polygonforclipping.shp out.shp in.shp . Users may have to install this first, but for the capstone they'll probably have to download Gdal (or maybe even QGis) anyways...

kdottiemo commented 7 years ago

@marwahaha: I could try replicating this. Spent yesterday trying a similar procedure with Python 3.4 and gave up during the waterfall of dependencies that I was encountering. If the importing process is frustrating for me, it's going to be deflating, at best, for the participants. Do you recommend foregoing Python 3.4 to try to make this capstone work with 2.7, or should we just use a regular Software Carpentry module as the fourth-quarter lesson at this point?

kdottiemo commented 7 years ago

Whoops, didn't mean to close.

marwahaha commented 7 years ago

Make an account

Make an account at NASA Earthdata. Follow this link for more info.

Download QGIS

We'll need this later. Download QGIS here.

Install pyModis and GDAL

We'll then install the Python package pyModis. pyModis is a python package to use when working with MODIS data. We can use it to download and convert data from EarthData. Try installing the package with pip, the Python package manager. Run the following:

conda install gdal
pip install pymodis

If the above doesn't work

Try installing the wheel by going here and downloading GDAL‑2.1.2‑cp27‑cp27m‑win_amd64.whl. Then, run pip install wheel and pip install <path/to/that/file>.
If that didn't work, on Windows and it complained about a missing vcvarsall.bat file, then download this compiler here from Microsoft and try again.

Get the scripts

We can use the pyModis to make our own Python scripts to download data, but someone else has already done that for us! In programming, there's a big culture of sharing code whenever you can. You should see scripts already in your path. Type

modis_downl

and press Tab. Does it tab-complete? Then you already have the scripts!

If the above doesn't work

You may need to download the scripts from the internet. Fortunately, this is all available on Github. Try cloning the repository with

git clone https://github.com/lucadelu/pyModis.git

This will make a folder called pyModis in your current directory. The scripts will be in the scripts/ folder of the pyModis.

Download MODIS data

Make a folder on your Desktop called MODIS, and cd into it. Make the 'MODIS' folder a git repository. Try doing all this on the command line.

We will download reflectance data of Japan over 8 days. More information about the product is available at this link.

Run the following, but replace the username and password with your EarthData account.

modis_download.py -U yourusername -P yourpassword -f 2012-06-01 -O -p MOD09A1.005 -t h28v05,h29v05,h28v04 .

Run ls in the current folder. Do you have several HDF files? If not, try repeating the command above.

About HDF

HDF stands for Hierarchical Data Format. It's a standard file type for data in many domains, including satellite data. It's great from an archivist's point of view --- you can store tons of metadata. It looks mostly like a set of excel tables, each with rows and columns.

If you're interested, you can view the HDF files by installing HDFView.

Make a Mosaic

Now, let's stitch together these HDF files to make a mosaic. This will lump together all the HDF files so we can work with each band.

modis_mosaic.py -s "1 0 1 1" -o mosaic -v listfileMOD09A1.005.txt

The -s flag determines which band to work with. In this case, we are working with bands 1, 3, and 4 (while ignorigng band 2).

This operation is quick. Run ls. Do you see files with a .vrt extension?

Good Practice

Try using ls with a wildcard to match all files that end in .vrt. Don't remember how to do that? Use ls *.vrt
Try using ls with arguments to print out more information about the file. What arguments do that? Use ls -al

Convert to GeoTIFF

Now, let's convert the VRT file to a GeoTIFF format. This is a standard image file (TIFF) with extra information encoded about its location. With GIS tools, you can use this as a layer in a map.

Let's convert one band, calling it final-01.tif:

modis_convert.py -v -s "( 1 )" -o final-01 -e 4326 A2012153_mosaic_sur_refl_b01.vrt

You should see a file called final-01.tif in the repository. You can open it up with a GeoTIFF viewer to see what the data shows, or even preview it with your Finder or File Explorer.

But there are multiple bands! Can you write a for loop to run through all the bands?

for filename in *.vrt
do
echo Converting $filename to $filename.tif ........................
modis_convert.py -v -s "( 1 )" -o $filename -e 4326 $filename
done

Removing the weird filename extension

Did your .tif files come out as .vrt.tif? This is because $filename corresponds to the entire input filename, including the extension. Hmm.... Sounds like we need to slice the name!

In Bash, you can add %.* to remove the last extension from a variable. So, if our files end in .vrt, the .vrt part will be removed.

Try running this:

for filename in *.vrt
do
echo Converting $filename to ${filename%.*}.tif ........................
modis_convert.py -v -s "( 1 )" -o ${filename%.*} -e 4326 $filename
done

Combining GeoTIFF files

We have a few different color bands of MODIS data, and want to combine the data into a single "color-composite" image. We can do that with another GDAL Python command.

gdal_merge.py -separate A2012153_mosaic_sur_refl_b01.vrt.tif A2012153_mosaic_sur_refl_b04.vrt.tif A2012153_mosaic_sur_refl_b03.vrt.tif -o stack.tif

This will take each band and stack them together in a new file called stack.tif.

Visualizing the GeoTIFF

For this part, you'll need QGIS. Download QGIS here.

Add a new layer (second button on the left sidebar), and select stack.tif. After some adjustement, you will see a true-color image of Japan. Pretty cool!

We will now clip our image to our study area extent, a common first step in our projects. We can use either the QGIS GUI or ... our command prompt! The cool thing about QGIS is that it dynamically writes the code as you change settings in the GUI tool.

Try this in the command line! Be sure to adjust your folder paths and filenames if they differ. gdalwarp -dstnodata -28762 -q -cutline C:/Users/Thinkpad/MODIS/studyarea_capstone/convex_japan_study_area.shp -crop _to_cutline -tr 0.00386756700055 0.00386760855155 -of GTiff C:/Users/Thinkpad/MODIS/stack2.tif clippedJapanRGB2

About Parameters

Try right-clicking on the layer, and selecting "Properties". Click on the "Style" tab, and you'll be able to adjust parameters that affect how the image is displayed. Oftentimes, your data isn't bad --- it's just not visualized well.

Challenge 1

Try to write a Bash script to download MODIS data, convert to .vrt, convert to GeoTIFF, and combine the color bands together. This time, use data from June 1st of 2013.

Challenge 2

Adjust the earlier Bash script so that it takes one argument: The date to download. You should be able to run sh Challenge2.sh 2010-06-01 to download and convert data from June 1st, 2010.

Challenge 3

Try to write a new Bash script to repeat this process, but download for June 1st of each year from 2005 to 2015. (Can you use a for loop with the script from the previous challenge?)

Summary

Here, we've done some real programming on some real data. We used Python to download actual MODIS data and convert it from standard HDF formats into GeoTIFFs. We assembled the different color bands into a single image, and displayed it using QGIS. We put this procedure in a Bash script so we can always run it, just specifying the date to download. Using the fundamentals of programming, we wrote loops to repeat this process over many years. We are using programming for Geospatial work in a real way, in just two days of learning.

But you all have just scratched the surface. As you keep coding, you will get stuck, and you will succeed, and you will get tired, and you will feel powerful. If you feel like you're repeating a task over and over, you can probably automate it. Talk to your fellow participants as you keep programming, and stay in touch. People are willing to help. Google is your resource. You are in charge of your computer.

You're becoming a programming wizard.

Thanks for coming to Software Carpentry. Now go analyze some data!