[R-sig-Geo] how to crete a *.prj file from a Raster object

[Roger Bivand]

And even:

require(raster)
require(rgdal)
r <- raster(list(x = 1:nrow(volcano), y = 1:ncol(volcano), z =
volcano), crs = "+proj=laea +lat_0=-50 +datum=WGS84")
writeRaster(r, "file.asc")
showWKT(proj4string(r), file="file.prj") 
# which also shows the WKT version on the console

Roger

--
Roger Bivand
NHH Norwegian School of Economics
Helleveien 30, 5045 Bergen, Norway
Roger.Bivand at nhh.no


________________________________________
Fra: r-sig-geo-bounces at r-project.org [r-sig-geo-bounces at r-project.org] på vegne av Michael Sumner [mdsumner at gmail.com]
Sendt: 31. juli 2013 10:20
Til: Emanuele Cordano
Kopi: r-sig-geo at r-project.org
Emne: Re: [R-sig-Geo] how to crete a *.prj file from a Raster object

One way is to write a dummy shapefile, and rename the .prj to suit.

To do that you could:
require(raster)
require(rgdal)
r <- raster(list(x = 1:nrow(volcano), y = 1:ncol(volcano), z =
volcano), crs = "+proj=laea +lat_0=-50 +datum=WGS84")
writeRaster(r, "file.asc")
dummy <- as(extent(r), "SpatialPolygons")
proj4string(dummy) <- CRS(proj4string(r))
writeOGR(SpatialPolygonsDataFrame(dummy, data.frame(x = 1)), ".",
"file", "ESRI Shapefile")
unlink(c("file.shp", "file.shx", "file.dbf"))


Some programs (such as GDAL) do use these .prj files for ASCII rasters
(see below), which helps make them more useful - though it's a bit of
a minefield since GDAL allows them to have non-square pixels, but the
raster package does not let you write one (rgdal::writeGDAL does). It
depends on your target for whether this mix is helpful or not. I would
seriously consider using a better format if you can. GeoTIFF is the
best all-rounder in my experience, but there's always lots of gotchas.
The native raster format (another .grd) is probably the best for
general R/raster package use if you don't need to go to other systems.

So, long story but just be aware that none of this makes for an
obvious choice, it depends on your circumstances. I avoid text-format
rasters like the plague.

On my system:

gdalinfo file.asc
Driver: AAIGrid/Arc/Info ASCII Grid
Files: file.asc
       file.prj
Size is 87, 61
Coordinate System is:
PROJCS["Lambert_Azimuthal_Equal_Area",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",-50],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (0.500000000000000,61.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Corner Coordinates:
Upper Left  (   0.5000000,  61.5000000) (  0d 0' 0.03"E, 49d59'58.01"S)
Lower Left  (   0.5000000,   0.5000000) (  0d 0' 0.03"E, 49d59'59.98"S)
Upper Right (  87.5000000,  61.5000000) (  0d 0' 4.39"E, 49d59'58.01"S)
Lower Right (  87.5000000,   0.5000000) (  0d 0' 4.39"E, 49d59'59.98"S)
Center      (  44.0000000,  31.0000000) (  0d 0' 2.21"E, 49d59'59.00"S)
Band 1 Block=87x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

rgdal can write these GDAL non-square pixels:

r <-  raster(volcano)
## this will fail
writeRaster(r, "afile.asc")
##Error in .startAsciiWriting(x, filename, ...) :
 ## x has unequal horizontal and vertical resolutions. Such data
cannot be stored in arc-ascii format

## this works
writeGDAL(as(r, "SpatialGridDataFrame"), "afile.asc", driver = "AAIGrid")

The crux of the GDAL rectangular pixel support is here at the top of the .asc:
ncols        61
nrows        87
xllcorner    0.000000000000
yllcorner    0.000000000000
dx           0.016393442623
dy           0.011494252874

This is described here, and how to overwrite the default:
http://gdal.org/frmt_various.html#AAIGrid

Cheers, Mike.

On Tue, Jul 30, 2013 at 11:53 PM, Emanuele Cordano
<emanuele.cordano at gmail.com> wrote:
> Dear list,
>
> I have a RasterLayer object with CRS.
>
> I would like to write it in a *.asc ESRI ASCII file with writeRaster. I did
> it , I got the *.asc file but no related *.prj file was written.
> Is there a quick method to create directly the *.prj file from a
> RasterLayer object?
>
> Thank you very much
>
>
> kind regards
> Emanuele Cordano
>
>
> --
> Emanuele Cordano
> Environmental Engineer / Ingegnere per l' Ambiente e il territorio nr.
> 3587 (Albo A - Provincia di Trento)
> Fondazione Edmund Mach
> CRI - Reasearch and Innovation Center
> Agro-Eco-System Dynamic Platform
> Via Mach 1
> 38017 San Michele All'Adige (Italy)
> phone +39 0461 615397
> email: emanuele.cordano at iasma.it ,emanuele.cordano 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



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

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo