I'm trying the following procedure (I apologize for mistakes but it is the very first time here). I've the raster dtm.temp expressed as follows:
> dtm.temp
class : RasterLayer
dimensions : 668, 965, 644620 (nrow, ncol, ncell)
resolution : 64.9, 92.6 (x, y)
extent : 437230.1, 499858.6, 5138842, 5200699 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:32632 +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : elev
values : 362.8404, 3584.865 (min, max)
Using this extent I convert it in longlat in the following way
# I build a spatial dataframe with extent data from dtm.temp
df <- data.frame(ID = 1:2, X = c(dtm.temp@extent@xmin, dtm.temp@extent@xmax),
Y = c(dtm.temp@extent@ymin, dtm.temp@extent@ymax))
coordinates(df) <- c("X", "Y")
crs_text <- crs(dtm.temp, asText=TRUE) # extracting crs from dtm.temp
proj4string(df) <- CRS(crs_text)
ext.lonlat <- spTransform(df, CRS("+proj=longlat +datum=WGS84"))
ext.lonlat
> ext.lonlat
class : SpatialPointsDataFrame
features : 2
extent : 8.183449, 8.998142, 46.40024, 46.95982 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
variables : 1
names : ID
min values : 1
max values : 2
at this point I use the extent (expressed as longlat) to crop the following raster dtm.temp1 (always in longlat but with a different initial extent and higher resolution)
> dtm.temp1
class : RasterLayer
dimensions : 9956, 14656, 145915136 (nrow, ncol, ncell)
resolution : 0.0002777778, 0.0002777778 (x, y)
extent : 7.857917, 11.92903, 44.27042, 47.03597 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : elev
values : -18, 4543 (min, max)
then
dtm.temp1.crop <- crop(dtm.temp1, ext.lonlat)
> dtm.temp1.crop
class : RasterLayer
dimensions : 2015, 2933, 5909995 (nrow, ncol, ncell)
resolution : 0.0002777779, 0.0002777772 (x, y)
extent : 8.183473, 8.998195, 46.40014, 46.95986 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
and then I project again in UTM the raster cropped at the same extent (or at least, so I believed..)
# crs_text as defined above
dtm.temp1.crop.proj <- projectRaster(dtm.temp1.crop, crs=crs_text)
>dtm.temp1.crop.proj
class : RasterLayer
dimensions : 2033, 2968, 6033944 (nrow, ncol, ncell)
resolution : 21.2, 30.9 (x, y)
extent : 437083.4, 500005, 5138362, 5201182 (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:32632 +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : elev
values : 360, 3603.655 (min, max)
as you can see the two extent dtm.temp (the starting one) e the final dtm.temp1.crop.proj seems to be too much different.
extent dtm.temp : 437230.1, 499858.6, 5138842, 5200699
extent dtm.temp1.crop.proj: 437083.4, 500005, 5138362, 5201182
The error (in meters) is:
> 437230.1 - 437083.4
[1] 146.7
> 499858.6-500005
[1] -146.4
> 5138842 - 5138362
[1] 480
> 5200699-5201182
[1] -483
I'm sure I'm doing something wrong but I'm going crazy trying to understand what...
Someone is so patient to help me? Thanks in advance!
I add some more information trying to clearify on the basis of the comment of Mike. Using a different approach (easier but computationally time consuming) I follow this steps:
> dtm.temp.extent <- extent(dtm.temp) # is the extent fixed
> dtm.temp.extent
class : Extent
xmin : 437230.1
xmax : 499858.6
ymin : 5138842
ymax : 5200699
# then I project directly the second raster dtm.temp1 using the crs
# with projectRaster (very slow)
> dtm.temp1.proj <- projectRaster(dtm.temp1,crs=crs_text)
# then I crop to the fixed extent
> dtm.temp1.proj.crop <- crop(dtm.temp1.proj, dtm.temp.extent)
# this is what I obtain
> extent(dtm.temp)
class : Extent
xmin : 437230.1
xmax : 499858.6
ymin : 5138842
ymax : 5200699
> extent(dtm.temp1.proj.crop)
class : Extent
xmin : 437233.6
xmax : 499852
ymin : 5138857
ymax : 5200688
the error now appears to me very reasonable:
> 437230.1 - 437233.6
[1] -3.5
> 499858.6 - 499852
[1] 6.6
> 5138842 - 5138857
[1] -15
> 5200699 - 5200688
[1] 11
The reason of the first approach is only because I'm trying to speedup the code (very time consuming in the second approach).
(Edited) I add, it may be useful for someone, an image of the problem in the first workflow (the two points of the extent - LL and UR - far from raster even if on a different raster respect to the one above in the question) and the solution suggested by Mkennedy that works fine (using four points unprojecting all four corners of the raster, LL, UL, UR, LR and then taking the min/max of the 8 values). I was not catching the rotation that occurs when unprojecting to lat/lon using only the LL and UR coordinates. Taking the true min/max, the results are more closely (as in the 2nd workflow)
Position of the raster vertex following a different approach
