[R-sig-Geo] resampling MODIS-based raster to PRISM raster to obtain same extent

Tue Jan 21 18:04:15 CET 2014

I think it would make more sense to use gdalwarp here, as you would be
done in one step (projection, resolution, and extent). Robert

On Tue, Jan 21, 2014 at 8:40 AM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
> Julie:
>
> First off, what you are attempting to do is a bit complicated, so try
> to not be intimidated!
>
> Second: it might help to understand that HDF files aren't rasters,
> persay, they are scientific data "containers" that can CONTAIN raster
> files.  So, often the first step is to extract the raster subcomponent
> of the HDF file into something more usable.
>
> As Robert mentioned, I have a package designed to make this quite a
> bit easier.  First, install a GDAL release that supports HDF4/5 from
> http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries.  On Windows,
> this is Osgeo4w.  On Mac this is the Kyngesburye GDAL Frameworks.
>
> Next:
>
> install.packages("gdalUtils", repos="http://R-Forge.R-project.org")
> # Note that I'm suggesting the r-forge version which has some key
> fixes vs. the CRAN version
>
> Third, you now can run the utilities.  You want to figure out which
> subset to extract first using:
> library("gdalUtils")
> ?get_subdatasets
> get_subdatasets("pathto/my.hdf") # Look at the examples of this function
> # This will return a list of subdatasets.  Make note of the number of
> the one you want.
>
> # Now, you can extract that subdataset to any format you want via
> gdal_translate:
> ?gdal_translate
> my_hdf_to_tiff <-
> gdal_translate("pathto/my.hdf","my_modis_extracted.tif",sd_index=X,output_raster=TRUE)
> # Where X is the number of the subdataset you got from get_subdatasets()
>
> # Finally, to sync all the files to the same projection, you can use
> my other package "spatial.tools"
> install.packages("spatial.tools")
> library("spatial.tools")
> # And use spatial_sync_raster
> ?spatial_sync_raster
>
> This function takes two inputs, "unsynced" (the file you are going to
> reproject/extend/crop) and "reference" (the file that serves as the
> reference projection, resolution, and extent).  Check the parameters
> to understand which resampling it will use.
>
> This puts a bunch of steps together in one: reprojection, resampling,
> cropping, and expanding the raster.  The final product will be the
> unsynced file that will exactly match the projection, extent and pixel
> size of the reference.
>
> Hope this helps!
>
> --j
>
> On Tue, Jan 21, 2014 at 7:42 AM, Julie Lee-Yaw <julleeyaw at yahoo.ca> wrote:
>> Thank you Agus and Robert for the suggestions. Unfortunately, I'm a little intimidated by GRASS and rgdal. I've had no experience with GRASS and very little experience with rgdal. Are there any examples you can point me to for the specific exercise of reprojecting raster data? I'm an absolute beginner here.
>>
>> Furthermore, the trick with the MODIS data is that the original downloaded data is in hdf format (which the raster package can't read…not sure about rgdal). Given this issue, I'm wondering if my strategy should be to initially reproject and save the MODIS data as geoTiff using the MRT tool and then use rgdal or GRASS to finish the job (e.g. specify the original reference system and then reproject to the parameters of the PRISM data)? Here I'm not even sure how to find information on the native origin of the MODIS files…I don't see anything in the MRT documentation about what the origin gets set to (you are right Robert…there  isn't an option for setting this so it must be doing it automatically).
>>
>> The PRISM data files only come with the following info in the prj file:
>>
>> GEOGCS["NAD83",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
>>
>> Is the PRIMEM["Greenwich", 0.0] specifying the origin in this case? Is origin the same as reference coordinate?
>>
>> Finally, in case it is relevant, both the PRISM and the MODIS rasters (after directly projecting using the MRT tool) give the following as the coord reference (and have the same resolutions)…
>>
>> R +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
>>
>>
>> The extents on the other hand are different.
>>
>> Thanks again for the continued help! It is much appreciated!
>>
>> Julie
>>
>>
>>
>>
>>
>> On Tuesday, January 21, 2014 9:03:03 AM, Agustin Lobo <alobolistas at gmail.com> wrote:
>>
>> Another option would be doing it through grass commands.
>> Agus
>>
>> On Mon, Jan 20, 2014 at 10:01 PM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
>>> Julie,
>>>
>>> Yes, in that case the projection (coordinate reference system -- CRS)
>>> will be the same. That is sufficient for vector data (points, lines,
>>> polygons), but not for raster data. For raster data you also need to
>>> match the resolution and origin such that the cells are aligned. With
>>> origin I refer to the location nearest to (0, 0) that the edge of a
>>> particular raster could come if only the extent is changed (resolution
>>> is fixed). If origins are not the same, the rasters are not aligned,
>>> and you cannot directly compare values for matching cells (as cells do
>>> not match). Even if two rasters are aligned, you may still need to
>>> crop/expand one or both such that they get exactly the same extent
>>> (bounding box).  resample can come to the rescue here but it is
>>> preferable to avoid it. See the raster package vignette for more info.
>>>
>>> I am under the impression that the MRT does not allow you to set an
>>> origin. If that is true, then it is an inadequate tool for changing
>>> the projection of raster data and I would use GDAL instead. On linux
>>> you can have rgdal with HDF5 support, but not on windows. But on
>>> windows you can use command line GDAL (from FWtools)  perhaps via the
>>> new gdalUtils package.
>>>
>>> Robert
>>>
>>>
>>> On Thu, Jan 16, 2014 at 11:48 PM, Julie Lee-Yaw <julleeyaw at yahoo.ca> wrote:
>>>> Hi Robert,
>>>>
>>>> Thanks for the response. I have done a bit more reading but am still
>>>> struggling to understand: if I set the datum (NAD83) and the resolution of
>>>> the NDVI data to match that of the prism data is that not sufficient for
>>>> compatible projections? Does the datum not provide the origin? The extent is
>>>> what I am now trying to adjust (e.g. the NDVI data is for all of North
>>>> America, the PRISM data is for the US only). Or do you mean something else
>>>> by extent?
>>>>
>>>> The prism data has a prj file that reads:
>>>>
>>>> GEOGCS["NAD83",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
>>>>
>>>> Thanks again for the help!
>>>>
>>>> Julie
>>>>
>>>>
>>>>
>>>>
>>>> On Tuesday, January 14, 2014 8:11:01 PM, Robert J. Hijmans
>>>> <r.hijmans at gmail.com> wrote:
>>>> Julie,
>>>> You raise an important question that is often overlooked. While it is
>>>> possible to use resample as you suggest; you would want to avoid it
>>>> because it leads to loss of data quality (although in practice this is
>>>> often minimal and irrelevant).
>>>>
>>>> You state that you "project it to the same geographic coordinate
>>>> system as PRISM precip. data". But that is an incomplete statement for
>>>> raster data. What you need to do is to project the NDVI data to the
>>>> raster definition used by PRISM. That includes the coordinate system,
>>>> but also the origin (or extent) and resolution of that raster. I do
>>>> not think MRT supports setting these parameters; in which case you
>>>> should not use it. You can use GDAL instead. On Linux this can be done
>>>> with rgdal; on windows you can use FWTools instead, or perhaps the new
>>>> gdalUtils package?
>>>>
>>>> Robert
>>>>
>>>> On Mon, Jan 13, 2014 at 3:39 PM, Julie Lee-Yaw <julleeyaw at yahoo.ca> wrote:
>>>>> Hi
>>>>>
>>>>> Using a combination of the scripts provided here (by Babak N.):
>>>>> http://r-gis.net/?q=ModisDownload and the MODIS reproject tool, I've finally
>>>>> managed to download NDVI data for North America and project it to the same
>>>>> geographic coordinate system as PRISM precip. data (e.g.
>>>>> http://www.prism.oregonstate.edu) for the same time period.
>>>>>
>>>>> I now want to stack these two rasters. My NDVI layer has a greater extent
>>>>> than the PRISM layer so I crop the former by the latter using:
>>>>>
>>>>> croppedNDVI<-crop(NDVI,prism)
>>>>>
>>>>>
>>>>> But when I look at the resulting raster, I see that the extent still
>>>>> doesn't line up. I think this is an issue with cells of the two rasters
>>>>> being "off" centre from each other. I can use the following to get them to
>>>>> align:
>>>>>
>>>>> adjustNDVI<-resample(croppedNDVI,prism)
>>>>>
>>>>>
>>>>> Now I can stack the "adjustNDVI" and "prism" layers as the extents match.
>>>>> But I am wondering whether this is valid? Why were the initial rasters
>>>>> misaligned in the first place given that I specified the same resolution and
>>>>> geographic coordinate system/datum when I processed the MODIS file? I'm
>>>>> grateful for any clarification!
>>>>>        [[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
>>
>>>>>
>>>>
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>         [[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