[R-sig-Geo] Problem Setting Projection of Rasters upon Import
Rainer M Krug
Rainer at krugs.de
Mon Jul 15 09:04:04 CEST 2013
"Robert J. Hijmans" <r.hijmans at gmail.com> writes:
> 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.
I would suggest to add a warning when trying to set a crs when it is
already supplied by the file. This would make the behavior clearer.
Rainer
>
> 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
<#secure method=pgpmime mode=sign>
--
Rainer M. Krug
email: RMKrug<at>gmail<dot>com
More information about the R-sig-Geo
mailing list