2

I'm trying to crop a raster dtm.temp1 (with higher resolution and in WGS84) at the extent of another dtm.temp (coarse resolution in UTM). Using different approch I reach different results

I try following two different approach as follows: 1.

 > 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
  1. 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 (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 of Mkennedy in comments that works fine (using four points).

Raster projected in UTM

  • 2
    When you warp the bounding box from one projection into another it also rotates. The corner coordinates of the source box are not at the corners of the envelope of that rotated box any more. It is wrong belief that the source envelope would round-trip back to original. – user30184 Sep 11 '17 at 15:54
  • Also, if both rasters have different resolution, is too difficult to have the same extent (if you make a corners rectification) – aldo_tapia Sep 11 '17 at 16:00
  • @aldo I don't expect to have the same extent obviously – EdGe Sep 11 '17 at 16:21
  • 2
    For the first workflow, try unprojecting all four corners of the raster, LL, UL, UR, LR and then take the min/max of the 8 values. You're not catching the rotation that occurs when unprojecting to lat/lon if you only do the LL and UR coordinates. If you take the true min/max, your results should more closely match the 2nd workflow. – mkennedy Sep 11 '17 at 18:19
  • 1
    On a larger scale, you even get the edges bended, with a new extent coordinate in the middle of the corners (especially if the original edge was in parallel with the axes). – AndreJ Sep 12 '17 at 5:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.