R/spatial_data_openstreetmap.R

In jae0/ecmgis: Package that interacts with stmv.

data_spatial_openstreetmap = function( what, p, redo=FALSE, fn=NULL, ... ) {

  # testing open street maps data ...

  p_add = list(...)
  if (length(p_add) > 0 ) p = c(p, p_add)
  i = which(duplicated(names(p), fromLast = TRUE ))
  if ( length(i) > 0 ) p = p[-i] # give any passed parameters a higher priority, overwriting pre-existing variable



  # -----------------------------

    if (what=="NS_coastline_openstreetmap") {
    # https://planet.openstreetmap.org/  https://planet.openstreetmap.org/historical-shapefiles/processed_p.tar.bz2 (projection not clear)
    # http://data.openstreetmapdata.com/coastlines-split-4326.zip (WGS84)

      fn = file.path( p$data.directory, "polygons", "NS_coastline_openstreetmap.rdata" )
      if ( !redo ) {
        if ( file.exists(fn) )  {
        load( fn )
        if (exists("internal.crs", p)) pg = spTransform( pg, sp::CRS(p$internal.crs) )
        return (pg)
      }
      }

      message( "Downloading .. you will need to manully unzip this one." )
      fn1 = "http://data.openstreetmapdata.com/land-polygons-complete-4326.zip"
      dir.local = file.path( p$data.directory, "polygons", "openstreetmap", "coastlines-split-4326" )
      dir.create( dir.local, recursive=TRUE, showWarnings=FALSE )
      fn.local = file.path( dir.local, basename( fn1 ) )
      shapefilename_osm = file.path( dir.local, "land_polygons.shp" )

      download.file( url=fn1, destfile=fn.local )
      unzip( fn.local, exdir=file.path( p$data.directory, "polygons", "openstreetmap" )  )


      # above is a messy polygon .. tidy up using the coarse bounds from stats canada
      ns = data_spatial( "NS_boundary_statscan_crude", p=p)

      out = raster::shapefile( shapefilename_osm )
      out = spTransform( out, sp::CRS(p$internal.crs) )
      out = gSimplify( out, tol=1)  # in meters

      pg = gIntersection( out, ns, drop_lower_td=TRUE, byid=TRUE )

      pg = as( pg, "SpatialPolygonsDataFrame" )
      pg$id = 1
      pg = try( gUnaryUnion(pg, id=pg@data$id), silent=TRUE  )

      # sum(gIsValid(pg, byid=TRUE)==FALSE) # check if any bad polys?
      # pg = gBuffer(pg, byid=TRUE, width=0)
      # plot(pg)
      pg = sp::spChFIDs( pg, "NovaScotia" ) #fix id's

      save(pg, file=fn, compress=TRUE)
      return(pg)
    }

  # ---------------------------------


      if (what=="NS_openstreetmaps_test") {
        # from open street maps data .. not fully working .. test only via osmdata
        fn = file.path( p$data.directory, "polygons", "NS_openstreetmaps.rdata" )
        if ( !redo ) {
          if ( file.exists(fn) )  {
            load( fn )
            if (exists("internal.crs", p)) pg = spTransform( pg, sp::CRS(p$internal.crs) )
            return (pg)
          }
        }

        install.packages("osmplotr")

        library("osmdata")
        library("osmplotr")
        library("maptools")

        library(sf)

        xlim=c(-66.5,-58)
        ylim=c(41,47.2)
        bb = c( xlim[1], ylim[1], xlim[2], ylim[2])


        o = getbb("Nova Scotia, Canada")

        Q <- opq(o)
        # Q <- add_osm_feature(Q, key='boundaries', value='administrative' )
        # Q <- add_osm_feature(Q, key='natural', value='lake' )
        Q <- add_osm_feature(Q, key = 'natural', value = 'water')

        Qxml <- osmdata_xml(Q)
        Qsp = osmdata_sp(Q, Qxml)

        bbox = get_bbox(o)
        map <- osm_basemap(bbox = bbox, bg = 'gray20')
        map <- add_osm_objects(map, y$osm_polygons, col = 'gray40')
        print_osm_map(map)


        # above is a messy polygon .. tidy up using the coarse bounds from stats canada
        ns = data_spatial( "NS_boundary_statscan_crude", p=p)

        out = raster::shapefile( shapefilename_osm )
        out = spTransform( out, sp::CRS(p$internal.crs) )
  #      out = gSimplify( out, tol=1)  # in meters

        pg = gIntersection( out, ns, drop_lower_td=TRUE, byid=TRUE )

        pg = as( pg, "SpatialPolygonsDataFrame" )
        pg$id = 1
        pg = try( gUnaryUnion(pg, id=pg@data$id), silent=TRUE  )

        # sum(gIsValid(pg, byid=TRUE)==FALSE) # check if any bad polys?
        # pg = gBuffer(pg, byid=TRUE, width=0)
        # plot(pg)
        pg = sp::spChFIDs( pg, "NovaScotia" ) #fix id's

        save(pg, file=fn, compress=TRUE)
        return(pg)
      }

  # --------------------------


      if (what=="NS_ocean_openstreetmapdata") {
        # from open street maps data
        fn = file.path( p$data.directory, "polygons", "NS_ocean_openstreetmapdata.rdata" )

        if ( !redo ) {

        if ( file.exists(fn) )  {
          load( fn )
          if (exists("internal.crs", p)) pg = spTransform( pg, sp::CRS(p$internal.crs) )
          return (pg)
        }
        }

        message( "Downloading .. you will need to manully unzip this one." )

        fn1 = "http://data.openstreetmapdata.com/water-polygons-split-4326.zip"
        dir.local = file.path( p$data.directory, "polygons", "openstreetmap", "water-polygons-split-4326" )
        dir.create( dir.local, recursive=TRUE, showWarnings=FALSE )
        fn.local = file.path( dir.local, basename( fn1 ) )
        shapefilename_osm = file.path( dir.local, "water_polygons.shp" )

        download.file( url=fn1, destfile=fn.local )
        unzip( fn.local, exdir=file.path( p$data.directory, "polygons", "openstreetmap" )  )

        # above is a messy polygon .. tidy up using the coarse bounds from stats canada
        ns = data_spatial( "NS_boundary_statscan_crude", p=p)

        out = raster::shapefile( shapefilename_osm )
        out = spTransform( out, sp::CRS(p$internal.crs) )
  #      out = gSimplify( out, tol=1)  # in meters

        pg = gIntersection( out, ns, drop_lower_td=TRUE, byid=TRUE )

        pg = as( pg, "SpatialPolygonsDataFrame" )
        pg$id = 1
        pg = try( gUnaryUnion(pg, id=pg@data$id), silent=TRUE  )

        # sum(gIsValid(pg, byid=TRUE)==FALSE) # check if any bad polys?
        # pg = gBuffer(pg, byid=TRUE, width=0)
        # plot(pg)
        pg = sp::spChFIDs( pg, "NovaScotia" ) #fix id's

        save(pg, file=fn, compress=TRUE)
        return(pg)
      }




}