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

[lbarqueira]

Hello,

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_1.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?

You can find the link to the example_1.xyz file, here

Thank you

[technocrat]

@lbarqueira

the z values are missing.

The z values are there, but they aren't a simple features geometry object. What to do about it depends on your intent. {sf} has limited abilities to encode 3-dimensional data.

library(data.table)
library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(terra)
#> Warning: package 'terra' was built under R version 4.3.1
#> terra 1.7.55
#> 
#> Attaching package: 'terra'
#> The following object is masked from 'package:data.table':
#> 
#>     shift

etd <- fread("/Users/ro/projects/demo/grist.csv")
setnames(etd,1:3,c("x","y","z"))
etd_df <- as.data.frame(etd)
etd_sf <- st_as_sf(x = etd_df, 
                       coords = c("x", "y"),
                       crs = 4326)

head(etd_sf)
#> Simple feature collection with 6 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -9.023816 ymin: 38.0144 xmax: -9.022335 ymax: 38.0144
#> Geodetic CRS:  WGS 84
#>         z                  geometry
#> 1 155.810 POINT (-9.023816 38.0144)
#> 2 155.555  POINT (-9.02352 38.0144)
#> 3 155.275 POINT (-9.023224 38.0144)
#> 4 154.990 POINT (-9.022928 38.0144)
#> 5 154.774 POINT (-9.022631 38.0144)
#> 6 154.536 POINT (-9.022335 38.0144)

^{Created on 2023-10-20 with reprex v2.0.2}