Calculate overlap between ATLANTIS-bgm file and grid from NPZD model for each grid/polygon combo.

[alketh]

Hi Michael,

I am currently working on coupling our ATLANTIS model with a well developed NPZD model for the NorthSea to test the phytoplankton, zooplankton and nutrient dynamics in our model. To do so I need to spatially overlap both models. I was wondering if you could give me some guidance.

The package "ecocat" (ECOHAM coupled to ATLANTIS) has two reference dataframes:

1. bgm_df --> positions of the Atlantis polygons
2. ecoham_layout --> midpoints of the ECOHAM grid cells

I figured the intersection is pretty straight forward with raster::intersect. Unfortunately, to do so the dataframes have to be converted to spatial objects. What is the best way to do this in your opinion?

@mdsumner

[mdsumner]

You can read the raw bgm file in with rbgm::bgmfile and then convert it so SpatialPolygonsDataFrame with rbgm::boxSpatial.

From a single data frame form you could use spbabel::sp though that requires that you name things specifically, i.e. "x", "y", "branch", "order", "island_". I can help with that and isolate the code if needed, it's not too hard to split the data.frame by part/polygon and recompose it as Spatial so it's not really necessary to use another package. (It's just something I'd prefer we were more systematic about).

If the grid cells are regular and interpretable as cell centres, then you can turn it into a raster with raster::rasterFromXYZ which has the obvious x, y, z column ordering. (If they aren't exactly regular you can work around that, but it also begs the question of why not just store the raster transform in place of the raw coordinates etc.) As a raster there is rasterToPolygons to build a SpatialPolygonsDataFrame (with yes every 4 corner-pixel expanded into 5-coordinate polygon rings).

I've found these topology overlays can still be a bit tricky, I've used intersect and cover successfully in the past - the key was having very careful data preparation up to the overlay steps. At any rate I'm happy to help find the best ways to do this if you have example code to get the data frames.

[mdsumner]

Oh, you can also use raster to get the approximate per-cell coverage within polygons, if you don't actually need the pieces - but I gather you want the pieces to work with down stream as well?

[alketh]

Thanks, that acutally helps a lot. The rbgm functions work like a charm. I was able ro convert the bgm file to a SpatialPolygonsDataFrame in a few seconds. Unfortunately, it seems that the grid-layout is not regular.

df <- raster::rasterFromXYZ(xyz = ecoham_layout) Error in raster::rasterFromXYZ(xyz = ecoham_layout) : x cell sizes are not regular

I suppose this is due to the curvature of the earth. The grids are actual regular via visual inspection when I plot them on a map. Please bare with me as I am new to the whole topic of rasters and so on. Is it possible to convert the long/lat information of the mid-point-grid to a SpatialPolygonsDataFrame? Or is it easier to calculate the intersections with simplified long/lat dataframes?

[mdsumner]

It depends, it might be just small numeric fuzz. Try raster::rasterFromXYZ(xyz = ecoham_layout, digits = 4) and maybe some smaller values.

Will have a closer look now

[alketh]

It does work with digits = 3. Thanks!

[mdsumner]

Hmm, have a look at plot(diff(RNetCDF::var.get.nc(nc_read, variable = "longitude"))) from https://github.com/alketh/ecocat/blob/master/data-raw/create-example-dfs.R

Notice too how raster doesn't complain when you go this way:

raster(system.file(package = "ecocat", "extdata/volume.nc"))

It's not exactly ideal, but it's pretty common. Glad you got it, if they really meant a non-regular grid of some sort that can be handled, but I think here it's totally fine. (It does happen, ROMS and CMIP are truly curvilinear but often grids are just stored with longlat when there's a clear affine transform from a projection).

Sorry, this is my rant trigger :)

[mdsumner]

But finally, doesn't this mean you can just raster() that netcdf and short-circuit a bit of coordinate munging? (I might have missed some stuff).

[alketh]

To be honest, I don't really care how to do it, I am mostly interested in the final result :)

Unfortunately,

atlantis_df <- rbgm::bgmfile(x = system.file(package = "ecocat", "extdata/NorthSea.bgm")) %>%
   rbgm::boxSpatial()

 ecoham_df <- raster::rasterFromXYZ(xyz = ecoham_layout, digits = 3) %>%
  raster::rasterToPolygons()

  test <- intersect_ecocat(atlantis_df, ecoham_df)

does result in NULL.

[alketh]

Sorry for the confusion... (I am getting a little bit tired, sry).

Of course I meant to say:

test <- raster::intersect(atlantis_df, ecoham_df)
Loading required namespace: rgeos 
Warning messages:
1: In raster::intersect(atlantis_df, ecoham_df) : non identical CRS
2: In raster::intersect(atlantis_df, ecoham_df) :
   polygons do not intersect

[mdsumner]

No worries, the problem is no projection metadata on your raster polygons. ;)

projection(ecoham_df) <- CRS("+init=epsg:4326")
intersect(atlantis_df, spTransform(ecoham_df, projection(atlantis_df)))

class       : SpatialPolygonsDataFrame 
features    : 2360 
extent      : 3339139, 4504114, 2946915, 4347280  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +a=6378137.0 +es=0.006694380022900787 +lon_0=10d00 +lat_0=52d00 +x_0=4321000.0 +y_0=3210000.0 
variables   : 12
names       : label, nconn,   botz,         area, vertmix, horizmix, insideX, insideY, .bx0, boundary, box_id, ecoham_id 
min values  :  Box0,     1, -100.0,   2848221108,       1,        1, 3408425, 2998658,    0,     TRUE,      0,       565 
max values  :  Box9,    17,    0.0, 111503422221,    

So you could transform the bgm polygons to 4326, but the result looks ok to me and I'd default to the native bgm projection in absence of any other rationale.

Note that this is a bit simpler, because the auto-projection pickup (though I haven't checked everything thoroughly ...):

atlantis_df <- rbgm::bgmfile(x = system.file(package = "ecocat", "extdata/NorthSea.bgm")) %>%
   rbgm::boxSpatial()
ecoham_df <-  raster(system.file(package = "ecocat", "extdata/volume.nc")) %>%
  raster::rasterToPolygons() %>% spTransform(projection(atlantis_df))

raster::intersect(atlantis_df, ecoham_df)

[mdsumner]

Good luck, it's def bed time here.

Happy to take any suggestions on rbgm etc, I really hate the names I gave the functions now "bgmfile" and "bgmfiles" was really stupid ...

[alketh]

Wow, so awesome! Thank you so much.

Yeah well, bgmfile definetely is a bit confusing... Such a shame that renaiming functions has so many repercussions.