Trying to read in pix4d outputs- specifically point clouds, orthophotos, dsm_xyz, and dtm.
I’m going to assume I need stars and sf.
library(sf)Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(stars)Loading required package: abindlibrary(ggplot2)
library(dplyr)
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(rayshader)# Load local functions
devtools::load_all()ℹ Loading galenRTypically the projdir would have a better name (likely similar/same as the project name)
projname <- '4m_upandback'
projdir <- file.path('data', 'pix4dout')dtmpath <- file.path(projdir, projname, '3_dsm_ortho', 'extras', 'dtm')read it in?
dtm <- read_stars(file.path(dtmpath, paste0(projname, '_dtm.tif')))That reads in, but isn’t plotting
useRaster = TRUE should be better, but it errors. Likely do to small rounding in the x and y coords (see https://search.r-project.org/CRAN/refmans/spatstat.geom/html/plot.im.html). I’m sure that could be cleaned up.
plot(dtm, useRaster = FALSE)downsample set to 1
Where is this stuff in space? is it in the right place?
# This doesn't work right- runs forever.
# library(tmap)
#
# tmap_mode("view")
#
# tm_shape(land) +
# tm_raster(dtm)orthopath <- file.path(projdir, projname, '3_dsm_ortho', '2_mosaic')read it in?
ortho <- read_stars(file.path(orthopath, paste0(projname, '_transparent_mosaic_group1.tif')))orthostars object with 3 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
Min. 1st Qu. Median Mean 3rd Qu. Max.
4m_upandback_transparent_mo... 0 0 0 3.11974 0 234
dimension(s):
from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]
band 1 4 NA NA NA NA 4 bands. What are they? I expected rgb, and the 4th is uniform. Seems to be alpha, based on some searching.
useRaster = TRUE works here. BUt not in quarto? that’s weird. Set to false so it renders.
plot(ortho, useRaster = FALSE)downsample set to 4
Color?
colmos <- st_rgb(ortho)colmosstars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
4m_upandback_transparent_mo...
Length:100000
Class :character
Mode :character
dimension(s):
from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]plot(colmos)downsample set to 1
Wow. That just worked. 4th band is likely alpha values.
The settings for this seem to be in the DSM and orthomosaic tab and resolution matches the ortho, as far as I can tell.
dsm_rpath <- file.path(projdir, projname, '3_dsm_ortho', '1_dsm')read it in?
dsm_r <- read_stars(file.path(dsm_rpath, paste0(projname, '_dsm.tif')))dsm_rstars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
Min. 1st Qu. Median Mean 3rd Qu.
4m_upandback_dsm.tif -24.73554 -24.70843 -24.68753 -24.68809 -24.67387
Max. NA's
4m_upandback_dsm.tif -24.62663 96962
dimension(s):
from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]Really not super clear why useRaster = TRUE isn’t working for these, but whatever.
plot(dsm_r, useRaster = FALSE)downsample set to 1
Those are funny z-values. What is the reference- clearly not elevation.
How many pixels? 1.0983136^{7}. 11 million is a lot.
This is an xyz point file.The settings are in the Additional Outputs tab, and we can set the grid spacing there. This was done with a spacing of 100cm. That implies these are centers, and this is a much coarser grid than in the dsm raster, which is at GSD scale.
dsm_xyzpath <- file.path(projdir, projname, '3_dsm_ortho', '1_dsm')dsm_xyz <- readr::read_csv(file.path(dsm_xyzpath, paste0(projname, '_dsm_1cm.xyz')),
col_names = c('x', 'y', 'z', 'r', 'g', 'b'))Rows: 1490037 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (6): x, y, z, r, g, b
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.assume the last 3 cols are rgb, use my function
dsm_xyz <- dsm_xyz %>%
mutate(hexcol = rgb2hex(r,g,b))Ignore spatial for a minute, can I just slam it in ggplot? using z as fill first. If we want the rgb, it’ll either be an orthomosiac or we’ll need to do something 3d with color overlay
dsm_xyzgg <- ggplot(dsm_xyz, aes(x = x, y = y, fill = z)) + geom_raster() + coord_equal()
dsm_xyzgg
That was really weirdly fast.
Can I use the rgb to see the orthomosiac?
dsm_xyzorthogg <- ggplot(dsm_xyz, aes(x = x, y = y, fill = hexcol)) + geom_raster() + coord_equal() + scale_fill_identity()
dsm_xyzorthogg
Would be cool to do a 3d map with color, maybe {rayshader}? Would really be nice to be able to rotate it etc.
That has color and height- might be good to feed to a ML algorithm to ID rocks using both sets of info and their relationships to each other.
This is supposedly a 100cm grid, but it has 1.5 million pixels, which is a lot less than the dsm raster at GSD resolution, but is nowhere near the difference I’d expect. It’s approximately 10x less, but GSD is < 1cm. And even at 1cm, I’d expect 10,000x fewer pixels. And those pictures are clearly not 1m pixels. What IS the spacing? Iterates fastest on x
dsm_xyz[2, 'x'] - dsm_xyz[1, 'x'] x
1 0.01# dsm_xyz[2, 'y'] - dsm_xyz[1, 'y']So, what is 0.01? That’s a cm, isn’t it?
Assuming the same crs as the raster dsm, the LENGTHUNIT is 1 meter. So, 0.01 would be 1cm, and that makes more sense.
st_crs(dsm_r)Coordinate Reference System:
User input: WGS 84 / UTM zone 55S
wkt:
PROJCRS["WGS 84 / UTM zone 55S",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 55S",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",147,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",10000000,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 144°E and 150°E, southern hemisphere between 80°S and equator, onshore and offshore. Australia. Papua New Guinea."],
BBOX[-80,144,0,150]],
ID["EPSG",32755]]It’s not exactly a translation from the raster, because the raster has all the NA off the edges, and a step of
diff(st_get_dimension_values(dsm_r, 'x')[1:2])[1] 0.00556So, the raster should have about 4x, but also includes the NAs.
sum(is.na(dsm_r[[1]]))[1] 5896746That’s not entirely it. Still about 266k more pixels in the rnn than expected once we drop the NA and adjust by pixel size.
#non-na in raster
rnn <- (2416*4546)-sum(is.na(dsm_r[[1]]))
#pixel difference- how many dsm_r pixels per dsm_xyz pixels?
pixdif <- (0.01/diff(st_get_dimension_values(dsm_r, 'x')[1:2]))^2
expected_rnn <- nrow(dsm_xyz)*pixdif
expected_rnn-rnn[1] -266385.2Can I make it geographic (sf)? Assume the same crs as the dsm_r
dsm_xyzsf <- st_as_sf(dsm_xyz, coords = c('x', 'y'), crs = st_crs(dsm_r))use the sf plot method? bad idea. takes forever. The points (next) are points, but it might make more sense to make these stars anyway- it is gridded.
plot(dsm_xyzsf[,'z'])
Try stars- first, confirm it is gridded- what’s the step?
Need to do a bit of cleanup- the x’s have negatives when they turnover y’s, and the y’s have a ton of zeros because they iterate slower.
rawx <- diff(pull(dsm_xyz[, 'x']))
cleanx <- rawx[rawx > 0]
# rounding error
difclean <- abs(cleanx-cleanx[1])
# huh- some are still out. by a lot. It seems systematic and strange
all(difclean < 1e-8)[1] FALSEwhich(difclean > 1e-8) [1] 191679 192281 192883 343775 344443 345109 375364 378772 378773
[10] 380145 381514 382200 443113 443822 483744 484458 485170 485881
[19] 486590 487296 527316 528041 530943 764429 765245 766058 766870
[28] 767680 768486 773330 774145 774961 775777 822605 823429 824254
[37] 825076 825896 877470 878310 879147 879982 884988 999784 1000608
[46] 1001431 1002252 1002256 1003070 1003881 1004687 1070957 1071726 1136369
[55] 1137116 1137859 1143776 1144521 1145270 1146022 1146781 1162799 1163554
[64] 1205080 1205808 1206534 1207256 1207976 1212992 1333223 1384324 1385469
[73] 1389993 1390552 1429954 1430434 1457519 1457948 1458373 1458795 1459216
[82] 1459636 1460053 1460467 1460877 1461283 1461685 1462082 1462475 1462863
[91] 1463247 1463626 1464003 1464376 1464744 1465107 1465466 1465821 1466172
[100] 1466520 1466861 1467199 1467534 1467865 1468192 1468513 1468831 1469145
[109] 1469455 1469760 1470062 1470360 1470655 1470946 1471233 1471516 1471796
[118] 1472073 1472346 1472616 1472881 1473144 1473402 1473656 1473905 1474150
[127] 1474392 1474629 1474862 1475092 1475316 1475537 1475754 1475966 1476174cleanx[which(difclean > 1e-8)] [1] 0.06 0.04 0.03 0.02 0.05 0.09 0.04 0.02 0.02 0.02 0.02 0.02 0.02 0.04 0.02
[16] 0.06 0.08 0.10 0.13 0.15 0.03 0.04 0.04 0.02 0.04 0.07 0.08 0.11 0.14 0.03
[31] 0.05 0.07 0.09 0.02 0.04 0.05 0.08 0.11 0.02 0.03 0.07 0.10 0.02 0.02 0.03
[46] 0.05 0.06 0.03 0.14 0.15 0.20 0.03 0.06 0.02 0.03 0.06 0.29 0.25 0.22 0.17
[61] 0.10 0.03 0.05 0.03 0.04 0.07 0.09 0.10 0.02 0.03 0.02 0.02 0.03 0.05 0.02
[76] 0.04 0.02 0.05 0.08 0.10 0.10 0.11 0.11 0.14 0.17 0.21 0.24 0.27 0.30 0.33
[91] 0.38 0.39 0.42 0.45 0.49 0.53 0.55 0.58 0.61 0.65 0.69 0.70 0.73 0.76 0.80
[106] 0.84 0.86 0.90 0.93 0.97 1.00 1.01 1.05 1.09 1.12 1.15 1.17 1.21 1.25 1.28
[121] 1.31 1.32 1.37 1.40 1.43 1.46 1.48 1.52 1.55 1.59 1.62 1.64 1.68 1.71 1.74If I plot that does it become clear? Because I took diffs and cut out the big steps, I can’t just attach it. Will need to modify in place.
dsm_xyz <- dsm_xyz |>
mutate(xdifs = x - lag(x, 1),
ydifs = y - lag(y, 1)) |>
mutate(xdifs = ifelse(ydifs == 0, xdifs, NA),
ydifs = ifelse(ydifs == 0, NA, ydifs))The steps are where there are weird things happening along angled edges.
dsm_xyz_x <- ggplot(dsm_xyz, aes(x = x, y = y, fill = xdifs)) + geom_raster() + coord_equal() + scale_fill_gradient(low = 'white', high = 'red')
dsm_xyz_x
Can we see it better with size?
dsm_xyz_xp <- ggplot(dsm_xyz, aes(x = x, y = y, color = xdifs, size = xdifs)) + geom_point() + coord_equal() + scale_color_gradient(low = 'white', high = 'red')
dsm_xyz_xpWarning: Removed 2370 rows containing missing values or values outside the scale range
(`geom_point()`).
Yeah, it’s also centered on corners.
Check y in the same way
dsm_xyz_yp <- ggplot(dsm_xyz, aes(x = x, y = y, color = ydifs, size = ydifs)) + geom_point() + coord_equal() + scale_color_gradient(low = 'white', high = 'red')
dsm_xyz_ypWarning: Removed 1487668 rows containing missing values or values outside the scale
range (`geom_point()`).
The ydifs are all 0.01, and we can see the rows stepover on the left edge.
I think it’s safe to ignore the few discontinuities.
Now, actually make it stars
dsm_xyzstars <- st_as_stars(dsm_xyz, crs = st_crs(dtm))plot(dsm_xyzstars['z'])
plot(dsm_xyzstars['hexcol'])
ggplot version- orthomosaic
dsm_xyzstarsgg <- ggplot() +
geom_stars(data = dsm_xyzstars, aes(fill = hexcol)) +
scale_fill_identity() +
theme(legend.position = 'none') +
coord_equal()
dsm_xyzstarsgg
ggplot version- elevation
dsm_xyzstarsgg_z <- ggplot() +
geom_stars(data = dsm_xyzstars, aes(fill = z)) +
colorspace::scale_fill_continuous_sequential(palette = 'Terrain 2') +
theme(legend.position = 'none') +
coord_equal()
dsm_xyzstarsgg_z
Is it faster to just use geom_raster? Would need to have a non-spatial df as we did above. This just errors when I try it.
dsm_xyzstarsgg_z_raster <- ggplot() +
geom_raster(data = dsm_xyzstars, aes(fill = z)) +
colorspace::scale_fill_continuous_sequential(palette = 'Terrain 2') +
theme(legend.position = 'none') +
coord_equal()
dsm_xyzstarsgg_z_rasterWhat about geom_sf on the point version? I plotted it above, but haven’t tried ggplot. This is more specifically what the point cloud is, below. Change everything to color from fill.
dsm_xyzsfgg_z <- ggplot() +
geom_sf(data = dsm_xyzsf, aes(color = z)) +
colorspace::scale_color_continuous_sequential(palette = 'Terrain 2') +
theme(legend.position = 'none') +
coord_sf()
dsm_xyzsfgg_z
Should work the same as dsm_xyz- this is also an xyz point file.
pcpath <- file.path(projdir, projname, '2_densification', 'point_cloud')pc <- readr::read_csv(file.path(pcpath,
paste0(projname,
'_group1_densified_point_cloud.xyz')),
col_names = c('x', 'y', 'z', 'r', 'g', 'b'))Rows: 645290 Columns: 6
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (6): x, y, z, r, g, b
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.pc <- pc %>%
mutate(hexcol = rgb2hex(r,g,b))Like the dsm_xyz, I should be able to use ggplot directly. But here, the points are the points in the cloud, not the evenly spaced points of the dsm_xyz. So geom_raster doesn’t work. Could use geom_tile I guess, but it’s more appropriate to just use geom_point.
pcgg <- ggplot(pc, aes(x = x, y = y, color = z)) +
geom_point() +
coord_equal() +
colorspace::scale_color_continuous_sequential(palette = 'Terrain 2')
pcgg
And the rgb for the ortho
pcorthogg <- ggplot(pc, aes(x = x, y = y, color = hexcol)) +
geom_point() +
coord_equal() +
scale_color_identity()
pcorthogg
We can make that an sf. Doesn’t make sense to make it a stars, since it’s points and not gridded. The dsm_xyz is technically points, but they’re gridded so can go either way.
pcsf <- st_as_sf(pc, coords = c('x', 'y'), crs = st_crs(dtm))Now we can plot
plot(pcsf['z'])
plot(pcsf['hexcol'])
And ggplot/geom_sf
pcsfgg_z <- ggplot() +
geom_sf(data = pcsf, aes(color = z)) +
colorspace::scale_color_continuous_sequential(palette = 'Terrain 2') +
theme(legend.position = 'none') +
coord_sf()
pcsfgg_z
Point size is really making this ugly, but that’s not super important right now.
pcsfgg_z <- ggplot() +
geom_sf(data = pcsf, aes(color = hexcol)) +
scale_color_identity() +
theme(legend.position = 'none') +
coord_sf()
pcsfgg_z
Ok, that all works. Need to clean it up and make a simple package.
Do the pixels line up between the various outputs? The dsm_xyz is points, each with x,y,z,r,g,b, so that will just match. But they aren’t technically a raster. The ortho, raster dsm, and dtm are rasters, but it’s not clear their pixels match. So, let’s do some investigation and possibly re-processing if necessary.
st_dimensions(dsm_r) from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]st_dimensions(ortho) from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]
band 1 4 NA NA NA NA Confirm they match
all(st_get_dimension_values(dsm_r, 'x') == st_get_dimension_values(ortho, 'x'))[1] TRUEall(st_get_dimension_values(dsm_r, 'y') == st_get_dimension_values(ortho, 'y'))[1] TRUECan we combine those into one stars? the catch is that the rgb vals aren’t on the same scale as the z. For the ortho, we have a pixel value for each of 4 bands (r,g,b, alpha). But z isn’t a band. So how would we do it? Make the bands attributes, probably. Or just call z a band?
To make the bands attributes and add z, this works. Will need to better define what’s needed later before deciding whether to then merge this back or what. And maybe we should be using hex (e.g. merging colmos and dsm_r).
orthosplit <- split(ortho, 'band') |>
setNames(c('r', 'g','b', 'alpha'))
dsm_r <- setNames(dsm_r, 'z')
orthodsm <- c(orthosplit, dsm_r)Merging into a dim- takes a while. and it’s unclear what to call this dimension or its values- ‘band’ isn’t right, and the values aren’t 0-255 for the z. I think don’t do this unless we’re really sure we want to.
orthodsm_dim <- orthodsm |>
merge()Does it actually make more sense to add a z dimension? Maybe? we can’t just merge the dsm_r, because then it has no attributes. Can we control the merge above? Maybe, but it’ll be a bit tricky. ignore until we have a clearer definition of what we need. The few obvious things I’ve tried have failed.
Whatever we do about the combination, the pixels match for the DSM raster and the orthophoto
These are unlikely to match-the dtm’s settings are done elsewhere and unless we were really particular to set them the same as the ortho, (which would require downgrading the ortho), they won’t match. AND, it’s unnecessary, since there’s a raster dsm that does match.
Looking at the from-to, they clearly don’t match, and looking at dimensions confirms it
st_dimensions(ortho) from to offset delta refsys point x/y
x 1 2416 263755 0.00556 WGS 84 / UTM zone 55S FALSE [x]
y 1 4546 5774223 -0.00556 WGS 84 / UTM zone 55S FALSE [y]
band 1 4 NA NA NA NA st_dimensions(dtm) from to offset delta refsys point x/y
x 1 1932 263755 0.00695 WGS 84 / UTM zone 55S FALSE [x]
y 1 3636 5774223 -0.00695 WGS 84 / UTM zone 55S FALSE [y]But, the DTM is some sort of smoothed thing, and is set in the Additional Outputs tab. We could make the ortho matched, but it’d be contrived (and pointless, given the existence of the raster dsm.
Again, I don’t really see the point of this turning the dsm_xyz into an image (raster) is already done in the raster dsm + ortho. We can merge if we want. Just not sure exactly how we’d want to present that.
Would be cool to do height with z and color with hexcol to actually map the stream in 3d with photo overlay. Pix4d does it, but would be nice to do here too.
plot_gg(dsm_xyzstarsgg, ggobj_height = dsm_xyzstarsgg_z, scale = 50)
render_snapshot()
Getting it to work directly isn’t happening for some reason, but gg is.
It can’t do the photo point matching on the moving surface. So the images are there but it can’t get parallax and so no z.
Have funs for xyz and funs for tifs. But call those within specific funs for each output (e.g. dsm_xyz and point cloud should have their own funs, and the dsm_xyz should allow returning a raster in addition to points, but pc shouldn’t). Similar for a standard set of plot returns