rterrarasterizing

Convert from an irregular cell sizes xyz file to raster. Missing Z value in the obtained raster


I have a xyz file that I want to convert to a raster using terra library. Although I have been able to convert to a raster the z values are missing.

The code is as following:

# Purpose: convert from an irregular cell sizes xyz file to raster
library(terra)
library(sf)
library(tidyverse)
library(data.table)

# loads a xyz file
etd <- fread("./exemplo.xyz")

class(etd) # "data.table" "data.frame"

#   V1       V2      V3
data.table::setnames(etd,1:3,c("x","y","z")) 

ras <- terra::rast(etd, type = "xyz", crs = 4326)
# Error: [raster,matrix(xyz)] x cell sizes are not regular
# Let see another approach.

etd_df <- as.data.frame(etd)

class(etd_df)

names(etd_df)
# [1] "x" "y" "z"

class(etd) # "data.table" "data.frame"


# crs = WGS84 (EPSG 4326)
# convert to spatial vector

etd_sf <- sf::st_as_sf(x = etd_df, 
                        coords = c("x", "y"),
                        crs = 4326)

plot(etd_sf)

# convert to SpatVector to use with terra
v_etd_sf <- terra::vect(etd_sf)

# convert sf to Spatraster
r1 <- terra::rast(v_etd_sf, ncols=700, nrows=700)

final_raster <- terra::rasterize(v_etd_sf, r1)

plot(final_raster)

range(final_raster$last) # min = 1, max = 1

# Questions:
# 1 - The Z has been lost;
# 2 - How to size, the values of ncols and nrows

Could you please help me find where is the problem?

File examplo.xyz can be find on the following link: exemplo.xyz

Thank you


Solution

  • You can interpolate into an xyz grid using interp. Here's a reproducible example:

    library(terra)
    library(interp)
    
    xyz <- read.table("original.xyz") 
    
    xyz$V3 <- as.numeric(sub(",", ".", xyz$V3))
    
    xyz <- interp::interp(x = xyz[[1]], y = xyz[[2]], z = xyz[[3]]) |> 
      interp::interp2xyz() |> 
      as.data.frame() |> 
      setNames(c("x", "y", "z"))
    
    ras <- rast(xyz, type = "xyz", crs = "WGS84")
    

    Resulting in

    ras
    #> class       : SpatRaster 
    #> dimensions  : 40, 40, 1  (nrow, ncol, nlyr)
    #> resolution  : 0.01738272, 0.01677493  (x, y)
    #> extent      : -10.00691, -9.311603, 37.99338, 38.66438  (xmin, xmax, ymin, ymax)
    #> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
    #> source(s)   : memory
    #> name        :        z 
    #> min value   : 2000.708 
    #> max value   : 3998.632 
    
    plot(ras)
    

    enter image description here