[R-sig-Geo] Problem Setting Projection of Rasters upon Import
Roger Bivand
Roger.Bivand at nhh.no
Fri Jul 12 10:00:01 CEST 2013
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]
>>>>>> 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
>>>>>>
>>>>>
>>>>>
>>>>
>>>> [[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
>>>
>>>
>>>
>>> --
>>> 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
>
--
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
More information about the R-sig-Geo
mailing list