"Robert J. Hijmans" <r.hijm...@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 <mtreg...@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, mtreg...@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 <mtreg...@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.biv...@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 <hodge...@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 [mtreg...@gmail.com] >>>>>> *Sent:* Thursday, July 11, 2013 8:21 PM >>>>>> *To:* Michael Sumner >>>>>> *Cc:* Hodgess, Erin; r-sig-geo@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 <mdsum...@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 <mtreg...@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 <mtreg...@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 <hodge...@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@r-project.**org<r-sig-geo-boun...@r-project.org>[ >>>>>>>>>> >>>>>>>>> r-sig-geo-bounces@r-project.**org <r-sig-geo-boun...@r-project.org>] >>>>>>> >>>>>>>> on behalf of Michael Treglia [mtreg...@gmail.com] >>>>>>>>>> Sent: Thursday, July 11, 2013 7:16 PM >>>>>>>>>> To: r-sig-geo@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@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@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: mdsum...@gmail.com >>>>>>> >>>>>>> >>>>>> >>>>>> >>>>> [[alternative HTML version deleted]] >>>>> >>>>> ______________________________**_________________ >>>>> R-sig-Geo mailing list >>>>> R-sig-Geo@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.biv...@nhh.no >>>> >>>> >>> >> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> R-sig-Geo mailing list >> R-sig-Geo@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 _______________________________________________ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo