2

I am working in flood forecasting and analysis with an intern 2D model which gives particular output of water level, not very good looking and not usable with meshers.

Today I use "SMS" a mesher from AQUAVEO which can interpolate data on the mesh. It can read the data that I extract from the model's results (CSV file that look like : X,Y,Z) for each time step I want to see :

> head(h)
    X       Y            Z
1 1005527 6280606 0.0029385281
2 1005537 6280614 0.0008433110
3 1005546 6280603 0.0009996744
4 1005552 6280620 0.0244638187
5 1005554 6280610 0.0008433110
6 1005556 6280594 0.0009992988

To resume I have n files for n time step of water level.

The mesher does not have a module that can assemble the time step and make an animation.

So end of the context! I wish to make a tool that can interpolate the data and make an animation of the results. For now I use R with the gstat package:

library(gstat)
library(raster)
library(sf)
library(sp)

h <- read.csv('Hmax.txt',sep = '',header = F,dec = '.',numerals = "no.loss")
crs <- "+init=EPSG:2154" #France
names(h) <- c("X","Y","Z")

bbox <- c(
  "xmin" = min(h$X),
  "xmax" = max(h$X),
  "ymin" = min(h$Y),
  "ymax" = max(h$Y)
)
sf_h <- st_as_sf(h, coords = c("X", "Y"), crs = CRS(crs))
nrows <- bbox["ymax"]-bbox["ymin"]
ncols <- bbox["xmax"]-bbox["xmin"]

h_rst <- raster(nrows = nrows, ncols = ncols,
                xmn = bbox["xmin"],ymx = bbox["ymax"],
                xmx = bbox["xmax"],ymn = bbox["ymin"], crs = CRS(crs))

fit_IDW <- gstat::gstat(
    formula = Z ~ 1,
    data = as(sf_h, "Spatial"),
    nmax = 16, nmin = 5,
    set = list(idp = i)
  )
  
  interp_IDW <- interpolate(h_rst,fit_IDW)

It works well ! But my model have 78 000 points! It takes 5 minutes by time step. (I know how to make the animation)

Do you know a way to make it quicker? Maybe with Python?

Jason Aller
  • 3,541
  • 28
  • 38
  • 38
Maxime auffret
  • 91
  • 5
  • Please make your code reproducible by providing the data in `dput` format. Please visit [How to make a great R reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example). – UseR10085 Oct 11 '21 at 10:41
  • Using `dput` with this type of data is very unwieldy and bad advice. You could instead build on the IDW example in `?raster::interpolate`. You do not say how many points you have, but if you have access to a HPC, you could speed things up by running many time steps in parallel. – Robert Hijmans Oct 11 '21 at 16:33
  • 1
    Thank you for your answer. I do say how many point I have : "It works well ! But my model have 78 000 points ! It takes 5 minutes by time step." . I could use HPC but the thing I don't understand is why R take so much time to interpolate (5 to 7 minutes) where Qgis/ArcGis or else will make it in few seconds. – Maxime auffret Oct 13 '21 at 09:18

0 Answers0