An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-sig-geo/attachments/20140430/7d4d7671/attachment.pl>
Problem merging rasters after conversion and reprojection of .hdf file
3 messages · Justin Michell, Albin Blaschka, Loïc Dutrieux
Hello! is the function resample() from the raster package what you need? HTH, Albin Am 30.04.2014 01:10, schrieb Justin Michell:
Hi Jonathan Thanks, very helpful! I have followed your suggestions (except I already had a loop going so I didn?t use batch_gdal_translate). But after this process when I I try to match these NDVI values (for Jan as shown below) to the coordinates of my dataset, I get mostly NA?s. I am not sure if this makes sense, although the reprojected map when I plot it looks fine I think. My coordinates in my spatial dataset represent 464 unique sample sites scattered across Kenya with a count of the number people examined and positive at each site, and downloaded from the malaria atlas project (http://www.map.ox.ac.uk/browse-resources/). Any thoughts would be appreciated! # Jan tiles (4 tiles cover Kenya)
for (i in 1 : 4){
+ # extract the NDVI band (sd_index=1) and converts to .tiff
+ gdal_translate(out.files[i], paste0("Jan", "_", i, ".tiff"), sd_index=1)
+ }
mosaic_rasters(gdalfile=c("Jan_1.tiff", "Jan_2.tiff", "Jan_3.tiff", "Jan_4.tiff"),dst_dataset="NDVI1_mosaic.tiff", separate=F, verbose=TRUE)
gdalwarp("NDVI1_mosaic.tiff","MODIS_NDVI/NDVIJan.tiff", s_srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",t_srs="+proj=longlat +datum=WGS84 +no_defs", srcnodata=-3000,dstnodata=-3000)
NDVI1 <- raster("MODIS_NDVI/NDVIJan.tiff")*0.0001
NDVI1 <- crop(NDVI1, extent(kenya))
NDVI1
class : RasterLayer dimensions : 1107, 953, 1054971 (nrow, ncol, ncell) resolution : 0.008397857, 0.008397857 (x, y) extent : 33.905, 41.90816, -4.671056, 4.625372 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : in memory names : NDVIJan values : -0.2, 0.9031 (min, max)
summary(NDVI1)
NDVIJan Min. -0.2000 1st Qu. 0.1813 Median 0.2892 3rd Qu. 0.4550 Max. 0.9031 NA's 652208.0000
extract(NDVI1, sitesPolyPoints at coords)
extract(NDVI1, sitesPolyPoints at coords) [1] 0.5094 NA 0.6557 NA NA 0.6575 0.6493 0.6469 NA 0.6607 NA NA NA 0.6814 NA NA 0.6620 NA NA NA NA 0.6749 0.6887 NA 0.6880 0.5959 NA [28] NA NA 0.6570 0.6212 NA NA NA NA NA 0.6647 NA NA 0.6296 NA NA 0.5831 NA NA NA 0.6147 NA 0.6167 NA 0.6258 0.6615 NA 0.6613 [55] NA 0.6626 NA NA 0.6406 NA 0.6451 0.6718 NA 0.6641 NA NA NA NA 0.6479 NA 0.6506 NA 0.6478 NA 0.6629 NA NA NA 0.6519 NA 0.6699 [82] NA NA NA NA 0.5996 NA 0.6694 NA 0.6855 NA 0.6625 NA NA NA NA NA NA NA 0.6389 NA NA NA NA NA NA NA NA [109] NA NA NA NA NA NA NA NA NA NA NA 0.6203 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [136] 0.6159 NA NA 0.6538 NA 0.6613 0.6731 NA NA NA 0.6781 NA 0.6090 NA NA 0.6054 NA NA NA 0.6810 NA NA NA NA NA NA NA [163] 0.6486 NA NA 0.6328 NA NA 0.6681 NA 0.6681 0.6766 NA NA 0.6910 NA NA NA NA NA NA NA NA NA 0.6706 NA NA NA NA [190] NA 0.6839 NA 0.6632 NA NA NA NA NA NA 0.6716 NA NA NA NA NA NA 0.6282 NA NA NA NA 0.6607 NA NA NA NA [217] NA NA NA NA NA 0.2105 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0.7046 NA 0.7367 NA 0.6968 NA NA [244] NA 0.5180 0.3300 0.7082 0.5610 0.3725 NA NA NA NA NA NA 0.2636 0.7032 0.6201 0.6996 NA NA 0.3071 0.3292 0.4628 NA NA NA NA NA NA [271] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [298] NA NA NA NA 0.7367 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [325] NA NA NA NA NA NA NA NA NA NA NA NA NA 0.7073 NA NA NA NA NA NA NA NA NA NA NA NA NA [352] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [379] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [406] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [433] NA NA NA NA NA NA NA NA NA 0.3652 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
head(sitesPolyPoints at coords)
coords.x1 coords.x2 [1,] 33.9940 0.1580 [2,] 34.1164 -0.0115 [3,] 34.1175 0.3770 [4,] 34.1277 -0.4222 [5,] 34.1318 -0.0035 [6,] 34.1350 0.4016 Cheers Justin On Apr 29, 2014, at 2:56 PM, Justin Michell <jwm302 at gmail.com> wrote:
Hi Jonathan Thanks for you help here. I get the following error when I tried the following as per you advice (and the example for mosaic_rasters). Do I not have the right GDALDriver?:
gdal_setInstallation()
valid_install <- !is.null(getOption("gdalUtils_gdalPath"))
valid_install
[1] TRUE
if(require(raster) && require(rgdal) && valid_install)
+ {
+ NDVI_1 <- mosaic_rasters(gdalfile=c(out.files[1],out.files[2], out.files[3],out.files[4], out.files[5], out.files[6]),dst_dataset="NDVI1_mosaic.hdf", separate=F, of=c("HDF4","HDF5"), verbose=TRUE)
+ gdalinfo("NDVI1_mosaic.hdf")
+ }
Checking gdal_installation...
GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalbuildvrt" -separate "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "MOD13A3.A2013001.h20v08.005.2013042085256.hdf" "MOD13A3.A2013001.h20v09.005.2013042085401.hdf" "MOD13A3.A2013001.h21v08.005.2013042085514.hdf" "MOD13A3.A2013001.h21v09.005.2013042085607.hdf" "MOD13A3.A2013001.h22v08.005.2013042085522.hdf" "MOD13A3.A2013001.h22v09.005.2013042085359.hdf"
Checking gdal_installation...
GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf"
Input file size is 3600, 24000
character(0)
attr(,"status")
[1] 1
Warning messages:
1: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf"' had status 1
2: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalinfo" "NDVI1_mosaic.hdf"' had status 1
ERROR 6: GDALDriver::Create() ... no create method implemented for this format.
ERROR 4: `NDVI1_mosaic.hdf' does not exist in the file system,
and is not recognised as a supported dataset name.
gdalinfo failed - unable to open 'NDVI1_mosaic.hdf?.
Cheers
Justin
On Apr 28, 2014, at 3:17 PM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Also note you can use ?batch_gdal_translate to do your HDF -> TIF conversions. --j On Mon, Apr 28, 2014 at 8:12 AM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Justin: Couple of comments: 1) In general, I find it easier to mosaic FIRST and then reproject SECOND -- when you reproject first, mosaic second, you can get weird issues particularly along the seams. 2) Since you are already using gdalUtils, have you tried ?mosaic_rasters (which will be faster than merge() in all likelihood, since it uses the base GDAL mosaicking capabilities)? Make sure separate=FALSE for this particular application. --j On Mon, Apr 28, 2014 at 7:18 AM, Justin Michell <jwm302 at gmail.com> wrote:
Dear Group I am dealing with the task of preprocessing and re-projecting MODIS gridded tiles in .hdf format from Sinusoidal projection, in R, to a geographic projection. I would like to align all my other climate layers (from worldclim.org) which are at 1km resolution to my NDVI layers. I am using the gdal_utils package in R, and only one subset I?m interested in (monthly NDVI). I am working on Mac OS and grabbed gdal from here: http://www.kyngchaos.com/software/frameworks . To do this I use the following code (previous list entries relating to gdalUtils have assisted me up till this point): # download all available monthly images for the year 2013 - MODIS PRODUCT (MOD13A3, Terra, Vegetation Indices, Tile, 1000m, Monthly) ModisDownload(x="MOD13A3",h=c(20,21,22),v=c(8,9),dates=c('2013.01.01','2013.12.31'),mosaic=F,proj=F) # extract the NDVI band (sd_index=1) and converts to .tiff for 1 tile for January gdal_translate("MOD13A3.A2013001.h20v08.005.2013042085256.hdf", "Jan_1.tiff", sd_index=1) # converts projection to Geog WGS84 gdalwarp(?Jan_1.tiff","MODIS_NDVI/NDVIJan_1.tiff", s_srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",t_srs="+proj=longlat +datum=WGS84 +no_defs",srcnodata=-3000,dstnodata=-3000) I have 6 tiles per month (covering Kenya). For each month I?d like to merge all 6 tiles using the raster package (to cover Kenya). For January I perform the following operation using the merge function in the raster package: NDVI1 <- (merge(raster("MODIS_NDVI/NDVIJan_1.tiff"),raster("MODIS_NDVI/NDVIJan_2.tiff"),raster("MODIS_NDVI/NDVIJan_3.tiff"), raster("MODIS_NDVI/NDVIJan_4.tiff"), raster("MODIS_NDVI/NDVIJan_5.tiff"), raster("MODIS_NDVI/NDVIJan_6.tiff"), tolerance=0.5))*0.0001 NDVI1 <- crop(NDVI1, extent(kenya)) Without including the tolerance argument I get an error: Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE, : different origin But after this merge a lot of NA?s are introduced and the plot doesn?t look right (as it did before for each tile). The resolution is not exactly the same. Here are a comparison of the layers (before merging NDVI tiles) and my SpatialPolygonsDataFrame called kenya, as well as my r sessionInfo() printed below:
n1 <- raster("MODIS_NDVI/NDVIJan_1.tiff")
n1
class : RasterLayer dimensions : 1219, 1275, 1554225 (nrow, ncol, ncell) resolution : 0.008205785, 0.008205785 (x, y) extent : 20, 30.46238, -0.002852273, 10 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : /Users/justinmichell/MODIS_NDVI/NDVIJan_1.tiff names : NDVIJan_1 values : -32768, 32767 (min, max)
meanTemp1
class : RasterLayer dimensions : 1115, 960, 1070400 (nrow, ncol, ncell) resolution : 0.008333333, 0.008333333 (x, y) extent : 33.90833, 41.90833, -4.666667, 4.625 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 data source : in memory names : layer values : -4.4, 30 (min, max)
kenya
class : SpatialPolygonsDataFrame features : 1 extent : 33.90722, 41.90517, -4.669618, 4.622499 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 variables : 11 names : FIPS, ISO2, ISO3, UN, NAME, AREA, POP2005, REGION, SUBREGION, LON, LAT min values : KE, KE, KEN, 404, Kenya, 56914, 35598952, 2, 14, 37.858, 0.53 max values : KE, KE, KEN, 404, Kenya, 56914, 35598952, 2, 14, 37.858, 0.53
sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gdalUtils_0.2.0 rgdal_0.8-14 RCurl_1.95-4.1 bitops_1.0-6 raster_2.2-12 sp_1.0-14
loaded via a namespace (and not attached):
[1] codetools_0.2-8 foreach_1.4.1 grid_3.0.2 iterators_1.0.6 lattice_0.20-24 scatterplot3d_0.3-34
Thanks very much,
Justin Michell
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 259 Computing Applications Building, MC-150 605 East Springfield Avenue Champaign, IL 61820-6371 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
-- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 259 Computing Applications Building, MC-150 605 East Springfield Avenue Champaign, IL 61820-6371 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
| Albin Blaschka, Mag.rer.nat. | Etrichstrasse 26, A-5020 Salzburg | * www.albinblaschka.info * | * www.researchgate.net/profile/Albin_Blaschka * | - It's hard to live in the mountains, hard but not hopeless!
Dear Justin,
gdalwarp is able to reproject/resample and mosaic all at once, and that
is probably the preferred option in that case (since there are no
overlapping areas).
See the simplified example below:
###
library(raster)
library(gdalUtils)
dir <- file.path(rasterOptions()$tmpdir, 'MODIS')
dir.create(dir, recursive=TRUE)
# Download two neighboring modis tiles
download.file('http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2000.04.22/MOD13A2.A2000113.h20v08.005.2006275210548.hdf',
destfile=file.path(dir, 'NDVI.h20v08.hdf'))
download.file('http://e4ftl01.cr.usgs.gov/MOLT/MOD13A2.005/2000.04.22/MOD13A2.A2000113.h21v08.005.2006275210549.hdf',
destfile=file.path(dir, 'NDVI.h20v09.hdf'))
# Get list of files downloaded
list <- list.files(dir, full.names=TRUE)
# Get list of sds1
sdsList <- sapply(X=list, FUN=function(x){get_subdatasets(x)[1]})
# Warp
gdalwarp(srcfile=sdsList, t_srs="+proj=longlat +datum=WGS84 +no_defs",
dstfile=file.path(dir, 'mosaic_latlong.tif'))
r <- raster(file.path(dir, 'mosaic_latlong.tif'))
plot(r)
Hope this helps
Best regards,
--
Lo?c Dutrieux
Laboratory of Geo-Information Science and Remote Sensing
Wageningen University
The Netherlands
On 14-04-30 01:10 AM, Justin Michell wrote:
Hi Jonathan Thanks, very helpful! I have followed your suggestions (except I already had a loop going so I didn?t use batch_gdal_translate). But after this process when I I try to match these NDVI values (for Jan as shown below) to the coordinates of my dataset, I get mostly NA?s. I am not sure if this makes sense, although the reprojected map when I plot it looks fine I think. My coordinates in my spatial dataset represent 464 unique sample sites scattered across Kenya with a count of the number people examined and positive at each site, and downloaded from the malaria atlas project (http://www.map.ox.ac.uk/browse-resources/). Any thoughts would be appreciated! # Jan tiles (4 tiles cover Kenya)
for (i in 1 : 4){
+ # extract the NDVI band (sd_index=1) and converts to .tiff
+ gdal_translate(out.files[i], paste0("Jan", "_", i, ".tiff"),
sd_index=1)
+ }
mosaic_rasters(gdalfile=c("Jan_1.tiff", "Jan_2.tiff", "Jan_3.tiff", "Jan_4.tiff"),dst_dataset="NDVI1_mosaic.tiff", separate=F, verbose=TRUE)
gdalwarp("NDVI1_mosaic.tiff","MODIS_NDVI/NDVIJan.tiff", s_srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",t_srs="+proj=longlat +datum=WGS84 +no_defs", srcnodata=-3000,dstnodata=-3000)
NDVI1 <- raster("MODIS_NDVI/NDVIJan.tiff")*0.0001
NDVI1 <- crop(NDVI1, extent(kenya))
NDVI1
class : RasterLayer dimensions : 1107, 953, 1054971 (nrow, ncol, ncell) resolution : 0.008397857, 0.008397857 (x, y) extent : 33.905, 41.90816, -4.671056, 4.625372 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : in memory names : NDVIJan values : -0.2, 0.9031 (min, max)
summary(NDVI1)
NDVIJan Min. -0.2000 1st Qu. 0.1813 Median 0.2892 3rd Qu. 0.4550 Max. 0.9031 NA's 652208.0000
extract(NDVI1, sitesPolyPoints at coords)
extract(NDVI1, sitesPolyPoints at coords) [1] 0.5094 NA 0.6557 NA NA 0.6575 0.6493 0.6469 NA 0.6607 NA NA NA 0.6814 NA NA 0.6620 NA NA NA NA 0.6749 0.6887 NA 0.6880 0.5959 NA [28] NA NA 0.6570 0.6212 NA NA NA NA NA 0.6647 NA NA 0.6296 NA NA 0.5831 NA NA NA 0.6147 NA 0.6167 NA 0.6258 0.6615 NA 0.6613 [55] NA 0.6626 NA NA 0.6406 NA 0.6451 0.6718 NA 0.6641 NA NA NA NA 0.6479 NA 0.6506 NA 0.6478 NA 0.6629 NA NA NA 0.6519 NA 0.6699 [82] NA NA NA NA 0.5996 NA 0.6694 NA 0.6855 NA 0.6625 NA NA NA NA NA NA NA 0.6389 NA NA NA NA NA NA NA NA [109] NA NA NA NA NA NA NA NA NA NA NA 0.6203 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [136] 0.6159 NA NA 0.6538 NA 0.6613 0.6731 NA NA NA 0.6781 NA 0.6090 NA NA 0.6054 NA NA NA 0.6810 NA NA NA NA NA NA NA [163] 0.6486 NA NA 0.6328 NA NA 0.6681 NA 0.6681 0.6766 NA NA 0.6910 NA NA NA NA NA NA NA NA NA 0.6706 NA NA NA NA [190] NA 0.6839 NA 0.6632 NA NA NA NA NA NA 0.6716 NA NA NA NA NA NA 0.6282 NA NA NA NA 0.6607 NA NA NA NA [217] NA NA NA NA NA 0.2105 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0.7046 NA 0.7367 NA 0.6968 NA NA [244] NA 0.5180 0.3300 0.7082 0.5610 0.3725 NA NA NA NA NA NA 0.2636 0.7032 0.6201 0.6996 NA NA 0.3071 0.3292 0.4628 NA NA NA NA NA NA [271] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [298] NA NA NA NA 0.7367 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [325] NA NA NA NA NA NA NA NA NA NA NA NA NA 0.7073 NA NA NA NA NA NA NA NA NA NA NA NA NA [352] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [379] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [406] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA [433] NA NA NA NA NA NA NA NA NA 0.3652 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
head(sitesPolyPoints at coords)
coords.x1 coords.x2 [1,] 33.9940 0.1580 [2,] 34.1164 -0.0115 [3,] 34.1175 0.3770 [4,] 34.1277 -0.4222 [5,] 34.1318 -0.0035 [6,] 34.1350 0.4016 Cheers Justin On Apr 29, 2014, at 2:56 PM, Justin Michell <jwm302 at gmail.com> wrote:
Hi Jonathan Thanks for you help here. I get the following error when I tried the following as per you advice (and the example for mosaic_rasters). Do I not have the right GDALDriver?:
gdal_setInstallation()
valid_install <- !is.null(getOption("gdalUtils_gdalPath"))
valid_install
[1] TRUE
if(require(raster) && require(rgdal) && valid_install)
+ {
+ NDVI_1 <- mosaic_rasters(gdalfile=c(out.files[1],out.files[2], out.files[3],out.files[4], out.files[5], out.files[6]),dst_dataset="NDVI1_mosaic.hdf", separate=F, of=c("HDF4","HDF5"), verbose=TRUE)
+ gdalinfo("NDVI1_mosaic.hdf")
+ }
Checking gdal_installation...
GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalbuildvrt" -separate "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "MOD13A3.A2013001.h20v08.005.2013042085256.hdf" "MOD13A3.A2013001.h20v09.005.2013042085401.hdf" "MOD13A3.A2013001.h21v08.005.2013042085514.hdf"
"MOD13A3.A2013001.h21v09.005.2013042085607.hdf" "MOD13A3.A2013001.h22v08.005.2013042085522.hdf" "MOD13A3.A2013001.h22v09.005.2013042085359.hdf"
Checking gdal_installation... GDAL command being used: "/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf" Input file size is 3600, 24000 character(0) attr(,"status") [1] 1 Warning messages: 1: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_translate" -of "HDF4" -of "HDF5" "/var/folders/t_/dby7tsnj0vgbfxp8hgclyrj80000gn/T//RtmpSKFW71/file17fd953983454.vrt" "NDVI1_mosaic.hdf"' had status 1 2: running command '"/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdalinfo" "NDVI1_mosaic.hdf"' had status 1 ERROR 6: GDALDriver::Create() ... no create method implemented for this format. ERROR 4: `NDVI1_mosaic.hdf' does not exist in the file system, and is not recognised as a supported dataset name. gdalinfo failed - unable to open 'NDVI1_mosaic.hdf?. Cheers Justin On Apr 28, 2014, at 3:17 PM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Also note you can use ?batch_gdal_translate to do your HDF -> TIF conversions. --j On Mon, Apr 28, 2014 at 8:12 AM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Justin: Couple of comments: 1) In general, I find it easier to mosaic FIRST and then reproject SECOND -- when you reproject first, mosaic second, you can get weird issues particularly along the seams. 2) Since you are already using gdalUtils, have you tried ?mosaic_rasters (which will be faster than merge() in all likelihood, since it uses the base GDAL mosaicking capabilities)? Make sure separate=FALSE for this particular application. --j On Mon, Apr 28, 2014 at 7:18 AM, Justin Michell <jwm302 at gmail.com> wrote:
Dear Group I am dealing with the task of preprocessing and re-projecting MODIS gridded tiles in .hdf format from Sinusoidal projection, in R, to a geographic projection. I would like to align all my other climate layers (from worldclim.org) which are at 1km resolution to my NDVI layers. I am using the gdal_utils package in R, and only
one subset I?m interested in (monthly NDVI). I am working on Mac OS and grabbed gdal from here: http://www.kyngchaos.com/software/frameworks .
To do this I use the following code (previous list entries relating to gdalUtils have assisted me up till this point):
# download all available monthly images for the year 2013 - MODIS PRODUCT (MOD13A3, Terra, Vegetation Indices, Tile, 1000m, Monthly)
ModisDownload(x="MOD13A3",h=c(20,21,22),v=c(8,9),dates=c('2013.01.01','2013.12.31'),mosaic=F,proj=F)
# extract the NDVI band (sd_index=1) and converts to .tiff for 1 tile for January
gdal_translate("MOD13A3.A2013001.h20v08.005.2013042085256.hdf", "Jan_1.tiff", sd_index=1)
# converts projection to Geog WGS84
gdalwarp(?Jan_1.tiff","MODIS_NDVI/NDVIJan_1.tiff", s_srs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs",t_srs="+proj=longlat +datum=WGS84 +no_defs",srcnodata=-3000,dstnodata=-3000)
I have 6 tiles per month (covering Kenya). For each month I?d like to merge all 6 tiles using the raster package (to cover Kenya). For January I perform the following operation using the merge function in the raster package:
NDVI1 <- (merge(raster("MODIS_NDVI/NDVIJan_1.tiff"),raster("MODIS_NDVI/NDVIJan_2.tiff"),raster("MODIS_NDVI/NDVIJan_3.tiff"), raster("MODIS_NDVI/NDVIJan_4.tiff"), raster("MODIS_NDVI/NDVIJan_5.tiff"), raster("MODIS_NDVI/NDVIJan_6.tiff"), tolerance=0.5))*0.0001
NDVI1 <- crop(NDVI1, extent(kenya))
Without including the tolerance argument I get an error:
Error in compareRaster(x, extent = FALSE, rowcol = FALSE, orig = TRUE, :
different origin
But after this merge a lot of NA?s are introduced and the plot doesn?t look right (as it did before for each tile).
The resolution is not exactly the same. Here are a comparison of the layers (before merging NDVI tiles) and my SpatialPolygonsDataFrame called kenya, as well as my r sessionInfo() printed below:
n1 <- raster("MODIS_NDVI/NDVIJan_1.tiff")
n1
class : RasterLayer dimensions : 1219, 1275, 1554225 (nrow, ncol, ncell) resolution : 0.008205785, 0.008205785 (x, y) extent : 20, 30.46238, -0.002852273, 10 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 data source : /Users/justinmichell/MODIS_NDVI/NDVIJan_1.tiff names : NDVIJan_1 values : -32768, 32767 (min, max)
meanTemp1
class : RasterLayer dimensions : 1115, 960, 1070400 (nrow, ncol, ncell) resolution : 0.008333333, 0.008333333 (x, y) extent : 33.90833, 41.90833, -4.666667, 4.625 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 data source : in memory names : layer values : -4.4, 30 (min, max)
kenya
class : SpatialPolygonsDataFrame features : 1 extent : 33.90722, 41.90517, -4.669618, 4.622499 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 variables : 11 names : FIPS, ISO2, ISO3, UN, NAME, AREA, POP2005, REGION, SUBREGION, LON, LAT min values : KE, KE, KEN, 404, Kenya, 56914, 35598952, 2, 14, 37.858, 0.53 max values : KE, KE, KEN, 404, Kenya, 56914, 35598952, 2, 14, 37.858, 0.53
sessionInfo()
R version 3.0.2 (2013-09-25)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gdalUtils_0.2.0 rgdal_0.8-14 RCurl_1.95-4.1 bitops_1.0-6 raster_2.2-12 sp_1.0-14
loaded via a namespace (and not attached):
[1] codetools_0.2-8 foreach_1.4.1 grid_3.0.2 iterators_1.0.6 lattice_0.20-24 scatterplot3d_0.3-34
Thanks very much,
Justin Michell
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 259 Computing Applications Building, MC-150 605 East Springfield Avenue Champaign, IL 61820-6371 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
-- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 259 Computing Applications Building, MC-150 605 East Springfield Avenue Champaign, IL 61820-6371 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo