[R-sig-Geo] Problem Setting Projection of Rasters upon Import
Robert J. Hijmans
r.hijmans at gmail.com
Sat Jul 13 01:29:21 CEST 2013
Michael,
I did not try it, but I think the below accomplishes the same. If at
all possible, you should avoid using ascii format files (that leads to
slow processing).
library(raster)
reference <- raster("D:/GIS/California/Hab Suitability
Files/10m_DEM/10m DEMasc/DEM_10m.asc")
projection(reference) <- "+proj=longlat +ellps=GRS80"
setwd("D:/GIS/PRISM/1981-2010/TMin")
s <- stack(list.files(pattern="*.asc", full.names=TRUE))
projection(s) <- "+proj=longlat +ellps=GRS80"
x <- crop(s, reference)
y <- projectRaster(x, crs="+proj=utm +zone=11 +datum=WGS84", res=800,
filename="TMinCrop.tif")
#if you need separate files you could use
writeRaster(y, .... , bylayer=TRUE)
As for setting the crs when creating a Raster object from file with
the raster function. I removed that because it is was some
inconsistent and misinterpreted. I will put it back in, allowing to
set the crs if it is unknown (as with esri ascii files), but not to
change it if the file supplies one.
Robert
On Fri, Jul 12, 2013 at 2:55 PM, Michael Treglia <mtreglia at gmail.com> wrote:
> Hi All,
>
> Thank you again for the suggestions and assistance. I used the work-around
> provided by Roger and Erin, and I've included the full function I was
> working on in case others find it useful.
>
> Basically, I've written this function to crop a bunch of ascii raster
> layers (.asc) in the working directory to the same extents as another
> raster layer (given that they are all in the same projection), then modify
> the projection and grain size of the cropped layers, as desired, and output
> them to new files in the working directory. In particular, I've been using
> this to reduce large raster grids from PRISM Climate Data to a focal area
> for which I need the data. The function can easily be modified to use
> different file types too...
>
> I'm all ears to improve the code... Depending on the size of your datasets
> it might take a few minutes to run. It does end up with a lot of data in
> memory for large files, so it would be great to hear any suggestions for
> dealing with that.
>
> #########################################
> #BatchCrop Function ###
> #by Mike Treglia, mtreglia at gmail.com ###
> ###Tested in R Version 3.0.0 (64-bit), using 'raster' version 2.1-48 and
> 'rgdal' version 0.8-10 ###
> ########################################
>
> #Function crops .asc raster files in working directory to extent of another
> layer (referred to here as 'reference' layer), converts to desired
> projection, and saves as new .asc files in the working directory. It is
> important that the original raster files and the reference layer are all in
> the same projection, though different pixel sizes are OK. Function can
> easily be modified to use other raster formats as well
>
> #Requires package 'raster'
>
> #Function Arguments:
> #'Reference' refers to name of the layer with the desired extent;
> 'OutName' represents the intended prefix for output files; 'OutPrj'
> represents the desired output projection; and 'OutRes' represents the
> desired Output Resolution
>
> BatchCrop<-function(Reference,OutName,OutPrj,OutRes){
> filenames <- list.files(pattern="*.asc", full.names=TRUE)
> #Extract list of file names from working directory
> library(raster) #Calls 'raster' library
>
> #Function 'f1' imports data listed in 'filenames' and assigns
> projection
> f1<-function(x,z) {
> y <- raster(x)
> projection(y) <- CRS(z)
> return(y)
> }
> import <- lapply(filenames,f1,projection(Reference))
>
> cropped <- lapply(import,crop,Reference) #Crop imported
> layers to reference layer, argument 'x'
>
> #Function 'f2' changes projectection of cropped layers
> f2<-function(x,y) {
> x<-projectRaster(x, crs=OutPrj, res=OutRes)
> return(x)
> }
> output <- lapply(cropped,f2,OutPrj)
>
> #Use a 'for' loop to iterate writeRaster function for all
> cropped layers
> for(i in (1:max(length(filenames)))){ #
> writeRaster(output[[i]],paste(deparse(substitute(OutName)),
> i), format='ascii')
> }
> }
> #############################################
> ###Example Code using function 'BatchCrop'###
> #############################################
>
> #Data layers to be cropped downloaded from:
> http://www.prism.oregonstate.edu/products/matrix.phtml?vartype=tmax&view=maps[testing
> was done using 1981-2010 monthly and annual normals; can use any
> .asc layer within the bounds of the PRISM data, with projection as lat/long
> and GRS80]
> #Set Working Directory where data to be cropped are stored
> setwd("D:/GIS/PRISM/1981-2010/TMin")
>
> #Import Reference Layer
> reference<-raster("D:/GIS/California/Hab Suitability Files/10m_DEM/10m DEM
> asc/DEM_10m.asc")
> #Set Projection for Reference Layer
> projection(reference) <- CRS("+proj=longlat +ellps=GRS80")
>
> #Run Function
> BatchCrop(Reference=reference,OutName=TMinCrop,OutPrj="+proj=utm +zone=11
> +datum=WGS84",OutRes=800)
>
>
>
> On Fri, Jul 12, 2013 at 11:39 AM, Michael Treglia <mtreglia at gmail.com>wrote:
>
>> Hi Roger,
>>
>> In R 2.13.0, it is Raster version 1.9-70; in R 3.0.0 it is Raster version
>> 2.1-48
>>
>> When running in R 3.0.0, upon running the 'raster' command, it loads the
>> rgdal package... whereas in R 2.13.0, I don't even have the rgdal package
>> installed. I guess in the previous version of Raster, it actually housed
>> gdal commands or something of the sort, whereas now it does not?
>>
>> Thanks for your help!
>> Mike
>>
>>
>> On Fri, Jul 12, 2013 at 1:00 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>>
>>> On Fri, 12 Jul 2013, Michael Treglia wrote:
>>>
>>> Thanks for trying Erin!
>>>> I appreciate knowing its not just me :-)
>>>>
>>>
>>> Are your raster package versions different? In current raster,
>>> R/rasterFromGDAL.R in line 129 sets crs="", and because in line 152 the
>>> assigned value in:
>>>
>>> projection(r) <- attr(gdalinfo, 'projection')
>>>
>>> is NA, you get what you see. My guess is that if your crs= argument is
>>> neither missing nor NA or the empty string, it should poverride
>>> attr(gdalinfo, 'projection').
>>>
>>> Providing the version of your older installation of raster will help
>>> isolate when and which change has led to this outcome, but you do have a
>>> workaround. I can't see the offending revision in SCM in R-forge.
>>>
>>> Roger
>>>
>>> -Mike
>>>>
>>>>
>>>> On Thu, Jul 11, 2013 at 6:40 PM, Hodgess, Erin <HodgessE at uhd.edu> wrote:
>>>>
>>>> Hi again.
>>>>>
>>>>> It looks like the crs argument no longer works when creating a raster
>>>>> from
>>>>> a file. I've tried all kinds of "variations on the theme" also, but no
>>>>> luck.
>>>>>
>>>>> Weird.
>>>>>
>>>>> Thanks
>>>>> Erin
>>>>>
>>>>> ------------------------------
>>>>> *From:* Michael Treglia [mtreglia at gmail.com]
>>>>> *Sent:* Thursday, July 11, 2013 8:21 PM
>>>>> *To:* Michael Sumner
>>>>> *Cc:* Hodgess, Erin; r-sig-geo at r-project.org
>>>>> *Subject:* Re: [R-sig-Geo] Problem Setting Projection of Rasters upon
>>>>>
>>>>> Import
>>>>>
>>>>> Hi,
>>>>> I apologize - I must have mistakenly copied/pasted only the last "t"
>>>>>
>>>>> Here is how it was run, and the results, in R 3.0.0 (64-bit):
>>>>>
>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="+proj=longlat
>>>>>> +ellps=GRS80")
>>>>>> test
>>>>>>
>>>>> class : RasterLayer
>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax, ymin,
>>>>> ymax)
>>>>> coord. ref. : NA
>>>>> data source : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>> names : us_ppt_1981_2010.02
>>>>> values : -2147483648, 2147483647 (min, max)
>>>>>
>>>>>
>>>>> and in R 2.13.0 (also 64-bit):
>>>>>
>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="+proj=longlat
>>>>>> +ellps=GRS80")
>>>>>> test
>>>>>>
>>>>> class : RasterLayer
>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax, ymin,
>>>>> ymax)
>>>>> coord. ref. : +proj=longlat +ellps=GRS80
>>>>> values : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Thu, Jul 11, 2013 at 6:13 PM, Michael Sumner <mdsumner at gmail.com
>>>>> >wrote:
>>>>>
>>>>> Please report on the actual code you ran, you've just sent two
>>>>>> messages, one with
>>>>>>
>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="+ellps=GRS80")
>>>>>>
>>>>>> this correctly cannot set the CRS and so it is NA
>>>>>>
>>>>>> and the other with
>>>>>>
>>>>>> t<-raster("us_ppt_1981_2010.**02.asc", crs="+proj=longlat
>>>>>> +ellps=GRS80")
>>>>>>
>>>>>>> test
>>>>>>>
>>>>>>
>>>>>> We have no idea what you see in "t" now. In short you have to have a
>>>>>> "+proj", not just a "+ellps" for this incantation. You can check with
>>>>>>
>>>>>> library(rgdal)
>>>>>> CRS("+ellps=GRS80")
>>>>>> Error in CRS("+ellps=GRS80") : projection not named
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> On Fri, Jul 12, 2013 at 11:10 AM, Michael Treglia <mtreglia at gmail.com>
>>>>>> wrote:
>>>>>>
>>>>>>> I just tried running it in an older version of R, 2.13.0 (64-bit) and
>>>>>>> it
>>>>>>> ran with no problem... the exact same code run from Notepad++
>>>>>>>
>>>>>>> Anybody know what has changed/if I need to adjust my for R 3.x.x?
>>>>>>>
>>>>>>> t<-raster("us_ppt_1981_2010.**02.asc", crs="+proj=longlat
>>>>>>> +ellps=GRS80")
>>>>>>>
>>>>>>>> test
>>>>>>>>
>>>>>>> class : RasterLayer
>>>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax,
>>>>>>> ymin,
>>>>>>> ymax)
>>>>>>> coord. ref. : +proj=longlat +ellps=GRS80
>>>>>>> values : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>>>>
>>>>>>> Thanks again!
>>>>>>> Mike
>>>>>>>
>>>>>>>
>>>>>>> On Thu, Jul 11, 2013 at 5:59 PM, Michael Treglia <mtreglia at gmail.com>
>>>>>>>
>>>>>> wrote:
>>>>>>
>>>>>>>
>>>>>>> Hi Erin,
>>>>>>>>
>>>>>>>> Thanks for the quick response - I did try it with the plus sign, and
>>>>>>>>
>>>>>>> same
>>>>>>
>>>>>>> result. (I also tried it with 'crs="+proj=longlat +ellps=GRS80"' with
>>>>>>>>
>>>>>>> the
>>>>>>
>>>>>>> same result too.)
>>>>>>>>
>>>>>>>>
>>>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="+ellps=GRS80")
>>>>>>>>> test
>>>>>>>>>
>>>>>>>> class : RasterLayer
>>>>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax,
>>>>>>>>
>>>>>>> ymin,
>>>>>>
>>>>>>> ymax)
>>>>>>>> coord. ref. : NA
>>>>>>>> data source : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>>>>> names : us_ppt_1981_2010.02
>>>>>>>> values : -2147483648, 2147483647 (min, max)
>>>>>>>>
>>>>>>>> I'm running R Version 3.0.0, 64-bit, in case that helps at all...
>>>>>>>> Cheers,
>>>>>>>> Mike
>>>>>>>>
>>>>>>>>
>>>>>>>> On Thu, Jul 11, 2013 at 5:50 PM, Hodgess, Erin <HodgessE at uhd.edu>
>>>>>>>>
>>>>>>> wrote:
>>>>>>
>>>>>>>
>>>>>>>> Hi Michael:
>>>>>>>>>
>>>>>>>>> Did you try this;
>>>>>>>>>
>>>>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="+ellps=GRS80")
>>>>>>>>>
>>>>>>>>> with the plus sign, please?
>>>>>>>>>
>>>>>>>>> That might do it.
>>>>>>>>>
>>>>>>>>> Thanks,
>>>>>>>>> Erin
>>>>>>>>>
>>>>>>>>> ______________________________**__________
>>>>>>>>> From: r-sig-geo-bounces at r-project.**org<r-sig-geo-bounces at r-project.org>[
>>>>>>>>>
>>>>>>>> r-sig-geo-bounces at r-project.**org <r-sig-geo-bounces at r-project.org>]
>>>>>>
>>>>>>> on behalf of Michael Treglia [mtreglia at gmail.com]
>>>>>>>>> Sent: Thursday, July 11, 2013 7:16 PM
>>>>>>>>> To: r-sig-geo at r-project.org
>>>>>>>>> Subject: [R-sig-Geo] Problem Setting Projection of Rasters upon
>>>>>>>>> Import
>>>>>>>>>
>>>>>>>>> Hi All,
>>>>>>>>>
>>>>>>>>> I am using he "raster" package, and in particular the "raster"
>>>>>>>>>
>>>>>>>> function. I
>>>>>>
>>>>>>> am dealing with a number of layers, and it would be easiest to set
>>>>>>>>>
>>>>>>>> the CRS
>>>>>>
>>>>>>> upon import. Looking through the documentation, it seems that it is
>>>>>>>>> doable,
>>>>>>>>> and I have been using the following code to get started, though it
>>>>>>>>> is
>>>>>>>>>
>>>>>>>> not
>>>>>>
>>>>>>> setting the CRS.
>>>>>>>>>
>>>>>>>>> Here are my code, and the results of the import:
>>>>>>>>>
>>>>>>>>>> test<-raster("us_ppt_1981_**2010.02.asc", crs="ellps=GRS80")
>>>>>>>>>> test
>>>>>>>>>>
>>>>>>>>> class : RasterLayer
>>>>>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax,
>>>>>>>>>
>>>>>>>> ymin,
>>>>>>
>>>>>>> ymax)
>>>>>>>>> *coord. ref. : NA *
>>>>>>>>> data source : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>>>>>> names : us_ppt_1981_2010.02
>>>>>>>>> values : -2147483648, 2147483647 (min, max)
>>>>>>>>>
>>>>>>>>> If I use the "projection" command afterwards, it sets the projection
>>>>>>>>> appropriately though.
>>>>>>>>>
>>>>>>>>>> projection(test)<-CRS("+proj=**longlat +ellps=GRS80")
>>>>>>>>>> test
>>>>>>>>>>
>>>>>>>>> class : RasterLayer
>>>>>>>>> dimensions : 3105, 7025, 21812625 (nrow, ncol, ncell)
>>>>>>>>> resolution : 0.008333333, 0.008333333 (x, y)
>>>>>>>>> extent : -125.0208, -66.47917, 24.0625, 49.9375 (xmin, xmax,
>>>>>>>>>
>>>>>>>> ymin,
>>>>>>
>>>>>>> ymax)
>>>>>>>>> *coord. ref. : +proj=longlat +ellps=GRS80 *
>>>>>>>>> data source : D:\GIS\PRISM\1981-2010\PPT\us_**ppt_1981_2010.02.asc
>>>>>>>>> names : us_ppt_1981_2010.02
>>>>>>>>> values : -2147483648, 2147483647 (min, max)
>>>>>>>>>
>>>>>>>>> Any suggestions of why the first option is not working? I have a
>>>>>>>>>
>>>>>>>> script in
>>>>>>
>>>>>>> which I import all files in a folder, and the simplest way to assign
>>>>>>>>>
>>>>>>>> the
>>>>>>
>>>>>>> projection would be the first option (I'm using an 'lapply' command).
>>>>>>>>> It would look something like this:
>>>>>>>>> import <- lapply(filenames ,raster, crs="....") #where "filenames"
>>>>>>>>> is
>>>>>>>>>
>>>>>>>> a
>>>>>>
>>>>>>> vector of filenames in the working directory.
>>>>>>>>>
>>>>>>>>> In case it helps, these are PRISM Climate Data, and the header of my
>>>>>>>>> raster
>>>>>>>>> file is this:
>>>>>>>>> ncols 7025
>>>>>>>>> nrows 3105
>>>>>>>>> xllcorner -125.020833333333329
>>>>>>>>> yllcorner 24.062500000000000
>>>>>>>>> cellsize 0.008333333333333
>>>>>>>>> NODATA_value -9999
>>>>>>>>>
>>>>>>>>> Thanks for any suggestions!
>>>>>>>>> Mike
>>>>>>>>>
>>>>>>>>> [[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<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<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>>>>>
>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> Michael Sumner
>>>>>> Hobart, Australia
>>>>>> e-mail: mdsumner at gmail.com
>>>>>>
>>>>>>
>>>>>
>>>>>
>>>> [[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<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>>
>>>>
>>> --
>>> Roger Bivand
>>> Department of Economics, NHH Norwegian School of Economics,
>>> Helleveien 30, N-5045 Bergen, Norway.
>>> voice: +47 55 95 93 55; fax +47 55 95 95 43
>>> e-mail: Roger.Bivand at nhh.no
>>>
>>>
>>
>
> [[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
More information about the R-sig-Geo
mailing list