Problem Setting Projection of Rasters upon Import

test

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

        [[alternative HTML version deleted]]

Thanks for trying Erin!
I appreciate knowing its not just me :-)

-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

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]]

        [[alternative HTML version deleted]]

 --
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com

	[[alternative HTML version deleted]]

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

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]]

        [[alternative HTML version deleted]]

 --
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com

        [[alternative HTML version deleted]]

--
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]]

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

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]]

        [[alternative HTML version deleted]]

 --
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com

        [[alternative HTML version deleted]]

--
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]]