Hi Robert,

Thanks a lot! That code worked like a charm - a lot shorter and more
elegant than my solution. Also, thanks for the explanation of why setting
CRS during the raster import was not working.

Have a great weekend!
Mike


On Fri, Jul 12, 2013 at 4:29 PM, Robert J. Hijmans <r.hijmans@gmail.com>wrote:

> 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@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@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@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@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@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@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 <mdsumner@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@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@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@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-bounces@r-project.org>[
> >>>>>>>>>
> >>>>>>>> r-sig-geo-bounces@r-project.**org <
> r-sig-geo-bounces@r-project.org>]
> >>>>>>
> >>>>>>> on behalf of Michael Treglia [mtreglia@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: mdsumner@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.Bivand@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
>

	[[alternative HTML version deleted]]

