Description Details Author(s) References See Also Examples
Process MODIS cloud mask product files to TIF, and then extract data
Package: | modiscloud |
Type: | Package |
Version: | 0.14 |
Date: | 2013-02-08 |
License: | GPL (>= 2) |
LazyLoad: | yes |
This package helps the user process downloaded MODIS cloud product HDF files to TIF, and then extract data. Specifically, MOD35_L2 cloud product files, and the associated MOD03 geolocation files (for MODIS-TERRA); and MYD35_L2 cloud product files, and the associated MYD03 geolocation files (for MODIS-AQUA).
The package will be most effective if the user installs
MRTSwath (MODIS Reprojection Tool for swath products;
https://lpdaac.usgs.gov/tools/modis_reprojection_tool_swath),
and adds the directory with the MRTSwath executable to
the default R PATH by editing ~/.Rprofile
.
Each MOD35_L2/MYD35_L2 file requires a corresponding MOD03/MYD03 geolocation file to be successfully processed with the MRTSwath tool.
MRTSwath is the MRT (MODIS Reprojection Tool) for the MODIS level 1 and level 2 products (cloud mask is level 2, I think).
A few example MODIS Cloud Product files, and derived
TIFs, are found in the data-only package
modiscdata
. These were too big to put in the main
package, according to CRAN repository policies
(http://cran.r-project.org/web/packages/policies.html).
Note: This code was developed for the following publication. Please cite if used: Goldsmith, Gregory; Matzke, Nicholas J.; Dawson, Todd (2013). "The incidence and implications of clouds for cloud forest plant water relations." Ecology Letters, 16(3), 207-314. DOI: http://dx.doi.org/10.1111/ele.12039
Nicholas J. Matzke matzke@berkeley.edu
https://lpdaac.usgs.gov/get_data/
NASA (2001). "MODLAND Product Filename Convention." <URL: http://landweb.nascom.nasa.gov/cgi-bin/QA_WWW/newPage.cgi?fileName=hdf_filename>.
Ackerman S, Frey R, Strabala K, Liu Y, Gumley L, Baum B and Menzel P (2010). "Discriminating clear-sky from cloud with MODIS algorithm theoretical basis document (MOD35)." MODIS Cloud Mask Team, Cooperative Institute for Meteorological Satellite Studies, University of Wisconsin - Madison. <URL: http://modis-atmos.gsfc.nasa.gov/_docs/MOD35_ATBD_Collection6.pdf>.
GoldsmithMatzkeDawson2013
check_for_matching_geolocation_files
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 | # Test function for checking roxygen2, roxygenize package documentation building
is.pseudoprime(13, 4)
# Some MODIS files are stored in this package's "extdata/" directory
# Here are some example MODIS files in modiscloud/extdata/
# Code excluded from CRAN check because it depends on modiscdata
## Not run:
library(devtools)
# The modiscdata (MODIS c=cloud data=data) package is too big for CRAN (60 MB); so it is available on github:
# https://github.com/nmatzke/modiscdata
# If we can't get install_github() to work, try install_url():
# install_github(repo="modiscdata", username="nnmatzke")
install_url(url="https://github.com/nmatzke/modiscdata/archive/master.zip")
library(modiscdata)
moddir = system.file("extdata/2002raw/", package="modiscdata")
# This directory actually has MYD files (from the MODIS-AQUA platform)
# (*will* work with the default files stored in modiscloud/extdata/2002raw/)
list.files(path=moddir, pattern="MYD")
# Check for matches (for MODIS-AQUA platform)
# (*will* work with the default files stored in modiscloud/extdata/2002raw/)
fns_df = check_for_matching_geolocation_files(moddir=moddir, modtxt="MYD35_L2", geoloctxt="MYD03", return_geoloc=FALSE, return_product=FALSE)
## End(Not run)
#######################################################
# Run MRTSwath tool "swath2grid"
#######################################################
# Source MODIS files (both data and geolocation)
# Code excluded from CRAN check because it depends on modiscdata
## Not run:
library(devtools)
# The modiscdata (MODIS c=cloud data=data) package is too big for CRAN (60 MB); so it is available on github:
# https://github.com/nmatzke/modiscdata
# If we can't get install_github() to work, try install_url():
# install_github(repo="modiscdata", username="nnmatzke")
# install_url(url="https://github.com/nmatzke/modiscdata/archive/master.zip")
library(modiscdata)
moddir = system.file("extdata/2002raw/", package="modiscdata")
# Get the matching data/geolocation file pairs
fns_df = check_for_matching_geolocation_files(moddir, modtxt="MYD35_L2", geoloctxt="MYD03")
fns_df
# Resulting TIF files go in this directory
tifsdir = getwd()
# Box to subset
ul_lat = 13
ul_lon = -87
lr_lat = 8
lr_lon = -82
for (i in 1:nrow(fns_df))
{
prmfn = write_MRTSwath_param_file(prmfn="tmpMRTparams.prm", tifsdir=tifsdir, modfn=fns_df$mod35_L2_fns[i], geoloc_fn=fns_df$mod03_fns[i], ul_lon=ul_lon, ul_lat=ul_lat, lr_lon=lr_lon, lr_lat=lr_lat)
print(scan(file=prmfn, what="character", sep="\n"))
run_swath2grid(mrtpath="swath2grid", prmfn="tmpMRTparams.prm", tifsdir=tifsdir, modfn=fns_df$mod35_L2_fns[i], geoloc_fn=fns_df$mod03_fns[i], ul_lon=ul_lon, ul_lat=ul_lat, lr_lon=lr_lon, lr_lat=lr_lat)
}
tiffns = list.files(tifsdir, pattern=".tif", full.names=TRUE)
tiffns
# For some unit testing etc., swath2grid may not be available. If so, use the default TIFs:
if (length(tiffns) == 0)
{
library(modiscdata)
tifsdir = system.file("extdata/2002tif/", package="modiscdata")
tiffns = list.files(tifsdir, pattern=".tif", full.names=TRUE)
}
#######################################################
# Load a TIF
#######################################################
library(rgdal) # for readGDAL
# numpixels in subset
xdim = 538
ydim = 538
# Read the grid and the grid metadata
coarsen_amount = 1
xdim_new = xdim / floor(coarsen_amount)
ydim_new = ydim / floor(coarsen_amount)
fn = tiffns[1]
grd = readGDAL(fn, output.dim=c(ydim_new, xdim_new))
grdproj = CRS(proj4string(grd))
grdproj
grdbbox = attr(grd, "bbox")
grdbbox
###########################
# Extract values from a particular pixel
###########################
# Greg's field site
greglat = 10.2971
greglon = -84.79282
grdr = raster(grd)
# Input the points x (longitude), then y (latitude)
point_to_sample = c(greglon, greglat)
xycoords = adf(matrix(data=point_to_sample, nrow=1, ncol=2))
names(xycoords) = c("x", "y")
xy = SpatialPoints(coords=xycoords, proj4string=grdproj)
#xy = spsample(x=grd, n=10, type="random")
pixelval = extract(grdr, xy)
# Have to convert to 8-bit binary string, and reverse to get the count correct
# (also reverse the 2-bit strings in the MODIS Cloud Mask table)
pixelval = rev(t(digitsBase(pixelval, base= 2, 8)))
print(pixelval)
## End(Not run)
#######################################################
# Load a TIF
#######################################################
# Code excluded from CRAN check because it depends on modiscdata
## Not run:
library(devtools)
# The modiscdata (MODIS c=cloud data=data) package is too big for CRAN (60 MB); so it is available on github:
# https://github.com/nmatzke/modiscdata
# If we can't get install_github() to work, try install_url():
# install_github(repo="modiscdata", username="nnmatzke")
# install_url(url="https://github.com/nmatzke/modiscdata/archive/master.zip")
library(modiscdata)
tifsdir = system.file("extdata/2002tif/", package="modiscdata")
tiffns = list.files(tifsdir, pattern=".tif", full.names=TRUE)
tiffns
library(rgdal) # for readGDAL
# numpixels in subset
xdim = 538
ydim = 538
# Read the grid and the grid metadata
coarsen_amount = 1
xdim_new = xdim / floor(coarsen_amount)
ydim_new = ydim / floor(coarsen_amount)
fn = tiffns[1]
grd = readGDAL(fn, output.dim=c(ydim_new, xdim_new))
grdproj = CRS(proj4string(grd))
grdproj
grdbbox = attr(grd, "bbox")
grdbbox
#######################################################
# Extract a particular bit for all the pixels in the grid
#######################################################
bitnum = 2
grdr_vals_bits = get_bitgrid(grd, bitnum)
length(grdr_vals_bits)
grdr_vals_bits[1:50]
#######################################################
# Extract a particular pair of bits for all the pixels in the grid
#######################################################
bitnum = 2
grdr_vals_bitstrings = get_bitgrid_2bits(grd, bitnum)
length(grdr_vals_bitstrings)
grdr_vals_bitstrings[1:50]
## End(Not run)
#######################################################
# Load some bit TIFs (TIFs with just the cloud indicators extracted)
# and add up the number of cloudy days, out of the total
# number of observation attempts
#######################################################
# Code excluded from CRAN check because it depends on modiscdata
## Not run:
library(devtools)
# The modiscdata (MODIS c=cloud data=data) package is too big for CRAN (60 MB); so it is available on github:
# https://github.com/nmatzke/modiscdata
# If we can't get install_github() to work, try install_url():
# install_github(repo="modiscdata", username="nnmatzke")
# install_url(url="https://github.com/nmatzke/modiscdata/archive/master.zip")
library(modiscdata)
tifsdir = system.file("extdata/2002bit/", package="modiscdata")
tiffns = list.files(tifsdir, pattern=".tif", full.names=TRUE)
tiffns
library(rgdal) # for readGDAL
# numpixels in subset
xdim = 538
ydim = 538
# Read the grid and the grid metadata
coarsen_amount = 1
xdim_new = xdim / floor(coarsen_amount)
ydim_new = ydim / floor(coarsen_amount)
sum_nums = NULL
for (j in 1:length(tiffns))
{
fn = tiffns[j]
grd = readGDAL(fn, output.dim=c(ydim_new, xdim_new))
grdr = raster(grd)
pointscount_on_SGDF_points = coordinates(grd)
grdr_vals = extract(grdr, pointscount_on_SGDF_points)
# Convert to 1/0 cloudy/not
data_grdr = grdr_vals
data_grdr[grdr_vals > 0] = 1
grdr_cloudy = grdr_vals
grdr_cloudy[grdr_vals < 4] = 0
grdr_cloudy[grdr_vals == 4] = 1
# Note: Don't run the double-commented lines unless you want to collapse different bit values.
# grdr_clear = grdr_vals
# grdr_clear[grdr_vals == 4] = 0
# grdr_clear[grdr_vals == 3] = 1
# grdr_clear[grdr_vals == 2] = 1
# grdr_clear[grdr_vals == 1] = 1
# grdr_clear[grdr_vals == 0] = 0
#
if (j == 1)
{
sum_cloudy = grdr_cloudy
#sum_not_cloudy = grdr_clear
sum_data = data_grdr
} else {
sum_cloudy = sum_cloudy + grdr_cloudy
sum_data = sum_data + data_grdr
}
}
# Calculate percentage cloudy
sum_nums = sum_cloudy / sum_data
grd_final = numslist_to_grd(numslist=sum_nums, grd=grd, ydim_new=ydim_new, xdim_new=xdim_new)
# Display the image (this is just the sum of a few images)
image(grd_final)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.