Couple of questions about buffers and MPAs in LSP

[mburgass]

Hi,

With LSP - I understand for the global it’s worked out as the amount of area that is protected compared to a reference point of 30%. What I was thinking about doing for this instead is to use the reference point as Ecologically or Biologically Significant Marine Areas and what % of these protected (reference point being 100% protection of these areas). I have a shape file of these areas and also would be using protectedplanet shapefile for the protected areas. I’m looking through the vector training module now. What I’m a little unsure of is how do I get buffers on my regions i.e. the 1km inshore buffer? Is this already done?

Also I’m a little bit confused as to the output – when I look on the LSP global data there is a km2 value for each year and each country all the way back to the 1800s. When you download the protectedplanet data do they have this through time? It looks like the values ping around a lot i.e. in 2003 it will be 626 and then same place in 2004 will be 292. Does this mean a huge chunk has been de-designated? Hope this makes sense. I think once I know how to do the buffer files I can start to work through this and it’ll probably become clearer.

[oharac]

Hi Mike,

Buffers

For the global analysis, we have buffer shapefiles for several offshore and inland distances, so it may be feasible to just use one of those, using something like rgdal::spTransform to convert it from its original projection to whatever projection you're using for your project.

However, it should be easy to generate a buffer (in case you want a custom distance, or if spTransform has problems projecting into a polar coordinate reference system). raster::buffer or rgeos::gBuffer can both take a polygon and extend it.

For offshore buffers, I've had success with this method: take the land shapefile and apply a buffer with a positive distance, and with byid = FALSE (for gBuffer; this creates a single feature with no region ID info. If you like raster::buffer I think it's dissolve = TRUE). Then intersect this feature (I like raster::intersect over rgeos::gIntersection here) with the EEZ polygons to chop it up into regions for analysis.

For inland buffers, you could do the same, from the other side - take the EEZ shapefiles, apply a buffer with byid = FALSE, and then intersect it with your land polygons (with their associated region IDs). If you use byid = TRUE, the buffers for each region will overlap; we don't have a pre-built way to cope with this yet.

Output for global LSP data

The old global method (e.g. v2015: https://github.com/OHI-Science/ohiprep/blob/master/globalprep/lsp/v2015/data/lsp_protarea_offshore3nm.csv) results in a dataframe with the amount of protected area added each year, rather than cumulative.

Starting in 2016 (v2016: https://github.com/OHI-Science/ohiprep/blob/master/globalprep/lsp/v2016/output/lsp_protected_inland1km.csv) the data_prep script calculates it as cumulative protected area, which seems more intuitive.

Whichever way you prefer, just make sure that your functions.R is interpreting the output correctly.

[mburgass]

Thanks Casey. Could you point me to the global buffers - I'll try this first?

Also I was wondering with creating MPA layers - what is the best way to do this? The protected planet dataset is very large. Is it a case of intersecting this with my regions? Thanks

[oharac]

Hi Mike,

The global buffers are on our NCEAS servers, so if you have access to Mazu you can find them. Here's the path to the region rasters once you're connected to Mazu (these are from the LSP script in globalprep/v2016/lsp_data_prep.Rmd):

• offshore 3nm: git-annex/globalprep/spatial/d2014/data/rgn_mol_raster_500m/rgn_offshore3nm_mol_500mcell.tif
• inland 1km: git-annex/globalprep/spatial/d2014/data/rgn_mol_raster_500m/rgn_inland1km_mol_500mcell.tif
• full eez: git-annex/globalprep/spatial/d2014/data/rgn_mol_raster_500m/rgn_eez_mol_500mcell.tif

The LSP data prep script has code for creating the MPA layers too. For global, we rasterize the MPA layers, which eliminates overlap (each 500 m cell has only one value, which is the first year under which it was protected).

• see rasterize_wdpa.py for a Python script we use (requires ArcPython) to rasterize the global WDPA dataset. Too tough for R to do in a reasonable amount of time. If you're good with Arc, go ahead and use this. Otherwise:
• After cropping down the global to a manageable size, you probably can get away with rasterizing in R; I'd suggest gdal_utils::gdal_rasterize() for speed and accuracy over raster::rasterize(), though the gdal function is trickier to use if you're not familiar with it.
  • I coded a more nicely-packaged version of gdal_rasterize, called gdal_rast2(). Feel free to copy it and adapt it as you see fit.
    • you tell it the path to the shapefile to rasterize (not a polygons object), hand it a base raster, give it a path to save the final raster, tell it which field in the shapefile to use for values.
    • If your projections don't match exactly it will give you an error; if you are sure that your projections match (the text lines may be different but refer to an identical coordinate reference system), you can override the proj4string check.
• Using polygons and intersecting with the MPA layers could result in some places double counting - e.g. SB Channel Islands would count as a National Park, marine protected areas, and National Marine Sanctuary. Plus intersecting, with many many polygons, will probably be painfully slow.

For a regional application, you're right, the global dataset is massive and hard to deal with. I would chop it down to size using something like:

• regional cropping - the raster::crop function works well to constrain a dataset to particular extents and works fine on polygons like the MPA dataset.
  • Create an extent object along the lines of arctic_ext <- extent(c(-180, 180, 60, 90) or whatever.
  • then crop the WDPA-MPA polygons object with something like arctic_mpa <- crop(global_wdpa_mpa, arctic_ext).
• filter by attributes such as country code, IUCN protected area classification, etc.
  • For global, we constrain to only "designated" protected areas (eliminating ones that are merely "proposed" etc)
  • additionally, we filter out some US polygons where the management plan (MANG_PLAN) is a non-MPA programmatic management plans - these are basically fisheries management areas, which don't really seem like "lasting special places"...
  • to see code for global, check out globalprep/v2016/lsp_prep_wdpa_poly.R

Note the coordinate reference systems here... the WDPA-MPA database is in WGS84, I think (degrees lat-long); but after filtering by attributes, we convert it to Mollweide for equal area calculations, then rasterize it. The region rasters are all in Mollweide projection so we can use the raster::crosstab function, and since all cells are the same area, we can just use cell counts, since we're just looking for proportion of protected area within the total area. You might not like Mollweide, but some type of equal-area projection would be easiest.

[mburgass]

Awesome thanks Casey. With regard to Mazu - I have access but when I login all I see is this:

Welcome to Mazu!

Login to RStudio Server at https://mazu.nceas.ucsb.edu/rstudio/.

Am I missing maybe access to the servers? Should I be able to see something else?

[mburgass]

Hi Casey,

I read in the TIF files to the rasters for the buffers, but I'm a bit unclear about how to use them. I tried to crop them / change the CRS but I think because they're read in as 'nativerasters' it doesn't like this command. Do you know how I get this in to something workable? Also it says it's 9.3gb - is doing anything with this in R just going to take forever?

Cheers,
Mike

[oharac]

Hi Mike,

I've never encountered nativeraster objects so not sure why they'd be showing up or what kind of trouble they might cause. Can you include your code for reading the file, any messages returned back, and info on the resulting raster? A quick peek at the documentation tells me that raster::raster() might return a nativeraster object if the native argument is TRUE (default is FALSE) or the rgdal package is missing, or for certain formats like NetCDF.

And yes, things will take a while... many raster package functions allow an argument progress = 'text' or similar so you get that running text bar to show how it's going.

cheers,

--Casey

[mburgass]

Hi Casey,

Back to this - I kind of abandoned this bit and moved on to something else but I've come back to it. I basically managed to do it all the way through for the whole EEZ as I had these spatial files. I then managed to do the 3nm buffer by expanding the land and overlapping with the EEZs like you said. I'm now on to the 1km buffer. Which I have taken from the global buffer file, cropped it down to the Arctic and reprojected. My issue is now how to separate it out in to regions - my issue being my land file is just one polygon. when I originally created my regions they were from offshore boundaries so onshore boundaries were never created. Is there any way to intersect the 1km inshore buffer with the EEZs to divide it up into regions?

Also I just wanted to check with my 3nm - I haven't quite finished it as when I built it, I actually went for what I think is 3km. I think we spoke about this on skype before but I can't quite remember. When using gbuffer as below - the width increases the buffer. If I set the width to 3000 - does that mean 3000m? I need to rejig this to make sure it is correct for 3nm. Thanks

crop_arc_land_buffer <- gBuffer(crop_arc_land, byid=FALSE, width=3000)

[oharac]

Hi Mike,

Re: units - If you check the projection/coordinate reference system of your polygons (use something like crop_arc_land@proj4string), you will probably see something in there that says '+units=m' which would tell you it's in meters (and 3 nm ends up being 5556 meters, I think).

Re: breaking inland buffer into regions - that's not straightforward to do based on the offshore EEZ regions. Your best bet is to get a land polygon with land borders (e.g. country borders, assuming that's how you're breaking it up) and intersect the inland 1 km buffer polygon with that (though if the coastlines don't match exactly you might end up with weirdness along there).

If you're good with GIS, you can also just draw big region blocks that extend farther inland and farther offshore than you need, but where the sides line up with the boundaries you want; when you intersect that with the buffer, the extra inland and offshore bits should disappear.

[mburgass]

Thanks Casey. I'm not good with GIS so think I will leave the 1km inland buffer for now. Not ideal, but too much other stuff to do. I'm encountering a weird error that is driving me crazy - When I try to use gbuffer for 5556m it causes R to crash - but only after about 20 minutes of processing or so. It worked fine with 3000m. I've tried deleting the repo, re-cloning, restarting R etc etc. I've tried it on my machine as well as mazu and the same issue happens. Could it be an issue with the further distance causing problems i.e. polygons overlapping?

[mburgass]

I'm going to give the 3nm raster from global a go - If I crop this down for the Arctic, do I just intersect it with EEZs to divide it up into regions? Thanks!

[mburgass]

Hi Casey- Sorry for all the messages :) I just tried the gbuffer with a value of 1000 - and it returned a 502 error but actually came through and worked. Seems to me like it's something to do with the bigger buffer area - but I'm not sure why this would make a difference? Unless it's timing out somehow?

[oharac]

hey Mike,

re: 5556 m buffer crashing, not sure why. It could be for the reason you described, or some other random issue. Is the "crash" a 502 error? because that may not be a crash - see below. If it is a legit crash, here are some options, in order of :

• instead of gBuffer(), try the buffer() function in the raster package, maybe with the dissolve = TRUE to get rid of the overlaps. Sometimes I would get one to work and not the other; they may have different methods for handling errors so the buffer may work for you where the gBuffer failed. Then you'd intersect that with your EEZ regions.
• use the 3nm buffer shapefile from global (using the wgs84 version if there's an option), create an extent object (from the raster package), something like arctic_ext <- extent(-180, 180, 60, 90) or whatever values, then crop (from raster package) the polygons to this extent (this will cut out most of it and make the later processing easier), probably do an spTransform to turn it from global CRS (likely WGS84) to your arctic CRS (some polar coordinate thing I suspect), and intersect that with your EEZ regions to divide it into your regions.
• use the 3nm raster from global, follow similar steps above - create an extent, crop it, then use projectRaster() instead of spTransform() to reproject the raster into your desired CRS.

The 502 error comes up when RStudio is expecting a response from R and doesn't get anything. This happens sometimes when you start running a long tedious process, and the RStudio window gets nervous that it hasn't received any news from the server for a while. Sometimes it's a crash, but often the process is still running in the background and will eventually report the results back to the RStudio window, and the 502 disappears. That's what it sounds like happened in your case.

Even once you get a proper buffer shapefile, rasterizing it can be tricky - I think I ended up resorting to arcPython to do the actual rasteribecause it was so challenging to get such large processes to work in R. Even just working at the arctic scale it still might be a big challenge. If you get stuck at this point, I can point you to an alternative raster function I wrote up (using gdalUtils::gdal_rasterize() instead of raster::rasterize) that might be better once you get to that point. Good luck!

Cheers,

--Casey