search

r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
564 stars 94 forks source link

write_mdim() stars object to view in QGIS as mesh layer #704

Open jvandens opened 3 months ago

jvandens commented 3 months ago

Consider the star object below, is there a way to save this using write_mdim() or other method that is compatible with viewing in QGIS as a Mesh Layer? https://docs.qgis.org/3.34/en/docs/user_manual/working_with_mesh/mesh.html#

For reference, this star object contains a polygon geometry stacked in z (depth) dimension and time dimension with 2 attributes (dataset groups in Mesh layer terminology)

image

star = structure(list(`F Coli` = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(geometry = 5L, 
z = 10L, times = 8L)), Entero = structure(c(1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1), .Dim = c(geometry = 5L, 
z = 10L, times = 8L))), dimensions = structure(list(geometry = structure(list(
    from = 1L, to = 5L, offset = NA_real_, delta = NA_real_, 
    refsys = structure(list(input = "NAD83 / UTM zone 18N", wkt = "PROJCRS[\"NAD83 / UTM zone 18N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    ID[\"EPSG\",26918]]"), class = "crs"), 
    point = FALSE, values = structure(list(structure(list(structure(c(498756.910095215, 
    510957.220092773, 522344.630126953, 512631.440124512, 498756.910095215, 
    4285731.00012207, 4304042.50012207, 4298020.00012207, 4279347.50012207, 
    4285731.00012207), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", 
    "sfg")), structure(list(structure(c(512631.440124512, 522344.630126953, 
    533818.380126953, 525469.630126953, 512631.440124512, 4279347.50012207, 
    4298020.00012207, 4292652.50012207, 4273846.50012207, 4279347.50012207
    ), .Dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), 
        structure(list(structure(c(525469.630126953, 533818.380126953, 
        545675.130126953, 538335.190124512, 525469.630126953, 
        4273846.50012207, 4292652.50012207, 4287713.50012207, 
        4268440.50012207, 4273846.50012207), .Dim = c(5L, 2L))), class = c("XY", 
        "POLYGON", "sfg")), structure(list(structure(c(538335.190124512, 
        545675.130126953, 558535.25012207, 552430.810119629, 
        538335.190124512, 4268440.50012207, 4287713.50012207, 
        4283077.50012207, 4262513.00012207, 4268440.50012207), .Dim = c(5L, 
        2L))), class = c("XY", "POLYGON", "sfg")), structure(list(
            structure(c(552430.810119629, 558535.25012207, 572934.75012207, 
            568734.440124512, 552430.810119629, 4262513.00012207, 
            4283077.50012207, 4279116.00012207, 4256515.00012207, 
            4262513.00012207), .Dim = c(5L, 2L))), class = c("XY", 
        "POLYGON", "sfg"))), class = c("sfc_POLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = 498756.910095215, 
    ymin = 4256515.00012207, xmax = 572934.75012207, ymax = 4304042.50012207
    ), class = "bbox"), crs = structure(list(input = "NAD83 / UTM zone 18N", 
        wkt = "PROJCRS[\"NAD83 / UTM zone 18N\",\n    BASEGEOGCRS[\"NAD83\",\n        DATUM[\"North American Datum 1983\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4269]],\n    CONVERSION[\"UTM zone 18N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-75,\n            ANGLEUNIT[\"Degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    ID[\"EPSG\",26918]]"), class = "crs"), n_empty = 0L)), class = "dimension"), 
    z = structure(list(from = 1L, to = 10L, offset = 1, delta = 1, 
        refsys = NA_character_, point = FALSE, values = NULL), class = "dimension"), 
    times = structure(list(from = 1L, to = 8L, offset = NA_real_, 
        delta = NA_real_, refsys = NA_character_, point = FALSE, 
        values = structure(list(start = c(0, 0.020833333954215, 
        0.0625, 0.10416666418314, 0.145833343267441, 0.1875, 
        0.22916667163372, 0.270833313465118), end = c(0.020833333954215, 
        0.0625, 0.10416666418314, 0.145833343267441, 0.1875, 
        0.22916667163372, 0.270833313465118, 0.3125)), class = "intervals")), class = "dimension")), raster = structure(list(
    affine = c(0, 0), dimensions = c(NA_character_, NA_character_
    ), curvilinear = FALSE, blocksizes = NULL), class = "stars_raster"), class = "dimensions"), class = "stars")
edzer commented 3 months ago

write_mdim() will write out the geometries to e.g. a netCDF file with discrete geometries (polygons); I don't know how the mesh stuff in qgis handles those, but it seems they could be quads. Any ideas or suggestions, @mdsumner ?

jvandens commented 3 months ago

reprex_model.nc.zip

Thanks @edzer ,if it helps @mdsumner , attached is the above star saved with write_mdim(). When loaded into QGIS, it seems to recognize some of the data that is in there, but is confused by the geometry and other dimensions (see first 2 screenshots)

For reference, the mesh are polygon quads, see last screenshot for what they should look like.

There should be 5 faces, 12 vertices repeated over 10 depths/elevations and 8 timesteps. With 2 datasets (entero and Fcoli) at each face.

image

image

image

jvandens commented 1 month ago

any thoughts on this?

edzer commented 1 month ago

Could you show (or sketch) a graph that shows what you would want to get?

jvandens commented 1 month ago

Hi @edzer I think the main issue is just a question of file format. I'm looking for a way that I can save a stars object into some format that can be read in as a "mesh layer" in QGIS (and from there, leverage all of the visualization tools in QGIS). This video has a really great overview of all those features, but it's missing the critical one of how to load the data in! From what I gather, QGIS mesh's should be able to read any of the mdal formats described here (but not all formats support all features). I think netCDF would do it but I'm not sure how it needs to be structured (and QGIS documentation doesn't seem to provide any examples and/or has broken links to them)

Here are two screenshots from that video that gives the idea. at 3:34: Mesh's support polygons with data described at faces and can be extruded over Depth/elevation and time dimensions (as my stars example in the OP) image

at 11:31: attached to each mesh can be any number of "datasets" which I think map to the attributes of the stars object. image

Thanks!