Thursday 30 August 2018

coordinate system - Reprojecting raster from lat/lon to UTM in R?


i have to turn it into a UTM in order to make the buffer functional.


wets<-readOGR(dsn=".",layer="shapefile")
r.raster <- raster()
extent(r.raster) <- extent(wets)
res(r.raster) <- 100

wets.r <- rasterize(wet,r.raster)
plot(wets.r)
wetsbuf<-buffer(wets.r,width=500)


During the buffer creation which is the last line of code, it gives this warning:


Warning message:  
In couldBeLonLat(x) :
raster has a longitude/latitude CRS, but coordinates do not match that

here's the info


  summary(wets.r)
layer
Min. 1

1st Qu. 1
Median 2
3rd Qu. 9
Max. 11
NA's 52629

summary(wets)

Object of class SpatialPolygonsDataFrame
Coordinates:

min max
x 683705 714088.8
y 4326266 4343768.0
Is projected: TRUE
proj4string :
[+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +datum=GGRS87
+units=m +no_defs +ellps=GRS80 +towgs84=-199.87,74.79,246.62]
Data attributes:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 2.5 5.0 5.0 7.5 10.0







wets.r

class : RasterLayer
dimensions : 175, 304, 53200 (nrow, ncol, ncell)

resolution : 100, 100 (x, y)
extent : 683705, 714105, 4326268, 4343768 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : layer
values : 1, 11 (min, max)
attributes :
ID FID
from: 1 0
to : 11 10


I have to change the prjection in order to be possible to do the buffer.



Answer



This is how you can reproject a raster in R using the raster package. In this example, the input geotiff was in a NAD83 geographic coordinate system and I reproject to a NAD 83 UTM 15 projected coordinate system. A good reference for Proj4 format projections, which are used by RGDAL, can be found at spatialreference.org.


library(raster)

# Create RasterLayer object
r <- raster('C:/temp/binary_nad83.tif')

# Define the Proj.4 spatial reference

# http://spatialreference.org/ref/epsg/26915/proj4/
sr <- "+proj=utm +zone=15 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"

# Project Raster
projected_raster <- projectRaster(r, crs = sr)

# Write the RasterLayer to disk (See datatype documentation for other formats)
writeRaster(projected_raster, filename="C:/temp/binary_utm15.tif", datatype='INT1U', overwrite=TRUE)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...