[GIS] How to get an ocean shapefile from a very detailed land area shapefile

clippolygonshapefile

I am looking for a possibility to "invert" a shapefile. What I want to do is to convert a shapefile with the continents and islands into a shapefile that has water polygons. The existing ocean shapefile from Natural Earth is too inaccurate.

As an active OSM contributor, I thought using the existing shapefiles from OSM would be fine for me. These are much more accurate but very buggy from a GIS perspective. However, mapserver can deal with it so I thought these should be fine. These shapefiles, some of them you can find here
http://metro.teczno.com/#coastline
are very large.

I have created a shapefile with one polygon with world bounds (-180 to 180, -90 to 90). Then I have tried to clip out the areas from my world bounds shapefile where land areas are in the other shapefile. I tried it with QGIS, vector tools -> difference function.
After 12 hours with 100% CPU usage, I stopped it. It always displayed a progress of 100%. RAM was not an issue.

I am looking for another possibility to clip out the land areas from a "world polygon", preferably Linux command line.

Best Answer

Rasters to the rescue! You can do this operation by converting your shapefile to a raster, doing the 'invert' operation and convert it back to vector

Raster -> Convert -> Rasterize. Output layer raster1
Open field calculator, enter expression raster1@1 = 0
Raster -> Convert -> Polygonize, Output layer ocean_polygons
Open Attribute table, search for DN = 0, toggle editing, delete selected features and save. What you will be left with will be polygons for the ocean areas.

enter image description here

Related Solutions

[GIS] Using a shapefile in QGIS makes it unusable in ArcGIS

QGIS also locks the shapefile. As long as it is opened in QGIS, you can not delete it from the file system (even if it is not in editing mode).

ARCGIS might not be pleased with that lock, and refuses to open the shapefile.

R – Removing Land Area from Marine Region Polygons

This might be a way.

I use this World map shapefile

# Read in new shapefile
world <- readOGR("~/test", layer="TM_WORLD_BORDERS-0.3")
# Order after $REGION in the 8th column
world_attri <- world@data[ order(-world@data[,8], world@data[,1]), ]
# Extract first line of each unique $REGION
world_attri <- world_attri[!duplicated(world_attri$REGION),]
# Subset the data to make the number of field equal to the Bahamian_shp
world_attri_subset <- subset(world_attri, select=c("NAME", "ISO3", "REGION"))
# Rename those fields to the same as the Bahamian_shp
colnames(world_attri_subset) <- c('name', 'id', 'mrgid')
# Reduce it by making a union based on $REGION
world_merged <- unionSpatialPolygons(world, world$REGION)
# Set the rownames to be equal to the mrgid field
rownames(world_attri_subset) <- world_attri_subset$mrgid
# Make a new SPDF
world_spdf <- SpatialPolygonsDataFrame(Sr=world_merged, data=world_attri_subset)

#### Following the above procedure ####

# Polygons seem to be plotted in alphabetical, so in order so make sure the bahamian ecoregion is plotted last we rename it
bahamian@data$name <- 'zBahamian'

# Insert the bahamian SPDF into the one for the world_spdf
new_spdf <- spRbind(bahamian_shp, world_spdf)
# ggplot requires a data frame, so use the fortify() function
mydf <- fortify(new_spdf, region = "name")

# Make a distinction between the underlying land shapes and the overlying ecoregion  polygon so that we can manually set the colours
mydf$filltype <- ifelse(mydf$id == 'zBahamian', "colour1", "colour2")

# Make the plot
ggplot(mydf, aes(x = long, y = lat, group = group)) +
  geom_polygon(colour = "black", size = 1, aes(fill = mydf$filltype)) +
  + coord_equal() +
  scale_fill_manual("Test", values = c(alpha("cyan", 0.4), "white"), labels = c("Bahamian", "World"))

### Calculate the corrected area (i.e. without land area) ###

diff <- gDifference(bahamian_shp, world_spdf, byid=FALSE, id=NULL)
gArea(SpatialPolygons(diff@polygons))
# 83.69948    ## Seems legit if one compares with the below area of the whole Bahamian

gArea(SpatialPolygons(bahamian_shp@polygons))
# 84.87498

(Area in square degrees since it's not projected in UTM)

enter image description here

This approach is probably more computationally demanding compared to the previous, but at least gaining efficiency in terms of reduced manual labour.

Related Question