Search code examples
rggplot2geospatialpolygonr-sf

fill polygon with points in R


I have a polygon which I would like to put some point on its perimeter and inside, with minimum distance as constraint.

here is my shape file :

    sh  = structure(list(farmId = "NO39", Name = "Sørlige Nordsjø I", 
    Country = "NO", Status = 99L, farmStatus = "Development Zone", 
    StatusComments = NA_character_, CapacityMWMin = 500, CapacityMWMax = 1500, 
    NoTurbinesMin = NA_integer_, NoTurbinesMax = NA_integer_, 
    Comments = "2013", 
    TurbineMWMin = NA_real_, TurbineMWMax = NA_real_, OtherNames = "Offshore wind - study area; Category A", 
    CountryName = "Norway", Lat = 57.4203547589749, Lon = 3.5335822891109, 
    IsEstimatedLocation = 0L, IsOnHold = 0L, OffshoreConsStartDate = structure(NA_real_, class = c("POSIXct", 
    "POSIXt"), tzone = ""), ProjFullCommDate = structure(NA_real_, class = c("POSIXct", 
    "POSIXt"), tzone = ""), Georegion = "Europe", Developers = "Norges vassdrags- og energidirektorat", 
    Owners = NA_character_, Operators = NA_character_, WaterDepthMinM = 58, 
    WaterDepthMaxM = 72, GISFILE = "EU_SCRE_WindFarms_PN_4C_230110_4258", 
    Shape_Length = 2.25250768130164, Shape_Area = 0.203134377892303, 
    Shape = structure(list(structure(list(list(structure(c(561491.358290986, 
    508749.11845824, 502715.786881769, 555417.318789534, 561491.358290986, 
    6359314.18214786, 6345317.40992537, 6369737.92848988, 6383257.96496235, 
    6359314.18214786), .Dim = c(5L, 2L)))), class = c("XY", "MULTIPOLYGON", 
    "sfg"))), n_empty = 0L, crs = structure(list(input = "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"Unknown based on WGS84 ellipsoid\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1],\n                ID[\"EPSG\",7030]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"UTM zone 31N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",3,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]],\n        ID[\"EPSG\",16031]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]", 
        wkt = "PROJCRS[\"unknown\",\n    BASEGEOGCRS[\"unknown\",\n        DATUM[\"Unknown based on WGS84 ellipsoid\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1],\n                ID[\"EPSG\",7030]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8901]]],\n    CONVERSION[\"UTM zone 31N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",3,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]],\n        ID[\"EPSG\",16031]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"), class = "crs"), class = c("sfc_MULTIPOLYGON", 
    "sfc"), precision = 0, bbox = structure(c(xmin = 502715.786881769, 
    ymin = 6345317.40992537, xmax = 561491.358290986, ymax = 6383257.96496235
    ), class = "bbox")), label = list(structure("Name: Sørlige Nordsjø I<br>Status: Development Zone<br>Capacity: 500[MW]<br>Owners: NA<br>", html = TRUE, class = c("html", 
    "character")))), sf_column = "Shape", agr = structure(c(WindfarmId = NA_integer_, 
Name = NA_integer_, Country = NA_integer_, Status = NA_integer_, 
WindfarmStatus = NA_integer_, StatusComments = NA_integer_, CapacityMWMin = NA_integer_, 
CapacityMWMax = NA_integer_, NoTurbinesMin = NA_integer_, NoTurbinesMax = NA_integer_, 
Comments = NA_integer_, TurbineMWMin = NA_integer_, TurbineMWMax = NA_integer_, 
OtherNames = NA_integer_, CountryName = NA_integer_, Lat = NA_integer_, 
Lon = NA_integer_, IsEstimatedLocation = NA_integer_, IsOnHold = NA_integer_, 
OffshoreConsStartDate = NA_integer_, ProjFullCommDate = NA_integer_, 
Georegion = NA_integer_, Developers = NA_integer_, Owners = NA_integer_, 
Operators = NA_integer_, WaterDepthMinM = NA_integer_, WaterDepthMaxM = NA_integer_, 
GISFILE = NA_integer_, Shape_Length = NA_integer_, Shape_Area = NA_integer_, 
label = NA_integer_), .Label = c("constant", "aggregate", "identity"
), class = "factor"), row.names = 1287L, class = c("sf", "data.frame"
))

And I could plot it using ggplot :

nc_l <- st_cast( sh , "MULTILINESTRING")

ggplot() +
    theme_void() +
    geom_sf(data = sh, fill = "grey") +
    geom_sf(data = nc_l, color = "red")

But when I'm tring to add the points I'm getting the CRS error ! the coorditanes are :

pts = structure(list(x = c(540775.687993986, 532383.520766688, 512939.996461768, 
521751.131317927, 545177.446866043, 513685.672432592, 523274.206243909, 
549810.734138201, 511952.141930947, 534994.719511726, 507718.99889467, 
546359.236299438, 537308.392091892, 509428.480302071, 523829.008167428, 
543606.738116335, 540456.878136299, 558469.563113825, 508799.852469142, 
551972.120242814, 547635.65103161, 541510.38374233, 517215.919804765, 
540521.856975332, 506263.078998375, 548914.157054223, 519569.151943714, 
547743.538814052, 532025.904750672, 514551.962500855, 510914.840154893, 
518107.052354577, 541002.347043822, 540081.536039312, 515329.281646527, 
554277.304895289, 545242.180912194, 542345.176843657, 521759.85060007, 
545181.936094893, 504533.562673026, 543532.513583723, 543621.478498706, 
529225.307147454, 520605.127557616, 525735.133265005, 508122.737869133, 
521682.429733122, 509955.311149239, 547794.732439182, 539651.510403828, 
558357.201812631, 544754.673356528, 554696.554826943, 535080.259461187, 
522320.18284748, 531092.594625451, 544112.304470398, 559613.707499225, 
536819.333643042, 533456.440618643, 561047.51031712, 559821.15251812, 
558594.60124343, 557367.85656924, 556140.91797012, 554913.78671035, 
553686.4622947, 552458.94478393, 551231.23425405, 550003.33017922, 
548775.23383893, 547546.94470733, 546318.46286027, 545089.78777157, 
543860.92072095, 542631.86118222, 541402.60923108, 540173.1649323, 
538943.52778143, 537713.69904724, 536483.67820311, 535253.46532457, 
534023.05988467, 532792.46316366, 531561.67463459, 530330.69437283, 
529099.52245372, 527868.15834991, 526636.60334192, 525404.85690239, 
524172.91909995, 522940.78942036, 521708.46913749, 520475.95772363, 
519243.25525382, 518010.36180308, 516777.27684329, 515544.00165552, 
514310.53571163, 513076.87908648, 511843.03125165, 510608.9934884, 
509374.76526573, 508492.58727485, 507972.78602581, 507453.50967794, 
506934.75779847, 506416.53176813, 505898.83115945, 505381.65614252, 
504865.00629269, 504348.88298289, 503833.28578741, 503318.21427692, 
502803.66982621, 503734.12515246, 504962.20241702, 506190.09474761, 
507417.80326881, 508645.32670649, 509872.66618499, 511099.82102979, 
512326.79056985, 513553.57592276, 514780.17641703, 516006.59137833, 
517232.82193024, 518458.8673983, 519684.72710819, 520910.40218312, 
522135.89134959, 523361.19573055, 524586.3146516, 525811.24743847, 
527035.99521375, 528260.55730306, 529484.93303218, 530709.12352331, 
531933.12810211, 533156.94609439, 534380.57862197, 535604.02441197, 
536827.28458596, 538050.35846965, 539273.24538891, 540495.94646493, 
541718.46102346, 542940.7883904, 544162.92968654, 545384.88423769, 
546606.65136976, 547828.23220319, 549049.62606378, 550270.83227752, 
551491.85196443, 552712.68385245, 553933.32906135, 555153.78691701, 
555832.1290514, 556361.64583031, 556891.67433984, 557422.21443943, 
557953.26592613, 558484.82865974, 559016.90303651, 559549.48771706, 
560082.58309663, 560616.18903572, 561150.30532924), y = c(6365554.63264462, 
6366094.06848788, 6352722.74197097, 6373206.96287922, 6364486.97954889, 
6349843.80155449, 6351633.30192643, 6372382.64955642, 6370706.02713503, 
6366631.00374243, 6365571.92404953, 6372105.23399408, 6366074.2155239, 
6370061.39157867, 6373624.33528036, 6364146.01192458, 6368459.81669427, 
6358706.91540676, 6361014.43382562, 6372167.4599021, 6381046.6079967, 
6369548.75553099, 6359274.09377177, 6358885.88908123, 6367143.40028134, 
6363328.44723995, 6352532.54030664, 6360915.23909248, 6369942.19520996, 
6364601.34093785, 6354541.65124067, 6369352.55366574, 6361244.67246172, 
6379102.41760098, 6353187.21342908, 6366030.12603715, 6373797.27609564, 
6378550.96244288, 6368486.02430351, 6376973.96349943, 6369502.33900287, 
6366652.91176626, 6375186.26077263, 6366649.38853267, 6355695.46125285, 
6353068.07546012, 6352965.38805807, 6364700.30686872, 6347236.9100541, 
6373250.16482407, 6354389.41572889, 6363055.9476361, 6371155.99216101, 
6359562.57826394, 6370816.83136007, 6354509.85888937, 6353852.73901094, 
6360315.62930042, 6366124.38073108, 6373785.09960235, 6369022.16867476, 
6359193.15701375, 6358859.05041212, 6358525.3127668, 6358191.94412723, 
6357858.94453459, 6357526.31516817, 6357194.05384269, 6356862.16172076, 
6356530.63885187, 6356199.48527816, 6355868.70106402, 6355538.28625142, 
6355208.24088991, 6354878.56502238, 6354549.25871182, 6354220.32200099, 
6353891.75493945, 6353563.55868997, 6353235.73107007, 6352908.27325416, 
6352581.18528595, 6352254.46721508, 6351928.11908619, 6351602.14095912, 
6351276.53287836, 6350951.29489359, 6350626.4270545, 6350301.9294067, 
6349977.80200829, 6349654.04490471, 6349330.65925891, 6349007.64289088, 
6348684.99697052, 6348362.72154404, 6348040.81666125, 6347719.28237196, 
6347398.11872352, 6347077.32577088, 6346756.90356124, 6346436.85214446, 
6346117.17156868, 6345797.86188744, 6345478.92426193, 6346349.69079351, 
6348443.00323905, 6350536.38618922, 6352629.84181068, 6354723.36893345, 
6356816.96638415, 6358910.63632997, 6361004.37648435, 6363098.18901566, 
6365192.07163692, 6367286.02651418, 6369380.05247518, 6369991.87262175, 
6370298.49973346, 6370605.49709366, 6370912.86465477, 6371220.60236559, 
6371528.71017931, 6371837.18804566, 6372146.03480044, 6372455.25262443, 
6372764.84035365, 6373074.7979366, 6373385.12532874, 6373695.8224788, 
6374006.88933481, 6374318.32585333, 6374630.13197934, 6374942.30767015, 
6375254.85287365, 6375567.76753706, 6375881.05161878, 6376194.70506624, 
6376508.72782612, 6376823.11985798, 6377137.88110872, 6377453.01152452, 
6377768.51106607, 6378084.37967451, 6378400.61731127, 6378717.22392241, 
6379034.19945329, 6379351.54386643, 6379669.25710741, 6379987.33912108, 
6380305.78987108, 6380624.60930249, 6380943.79735963, 6381263.3540073, 
6381583.27919004, 6381903.57285169, 6382224.23495816, 6382545.2654457, 
6382866.66428097, 6383188.43140768, 6381614.08159422, 6379517.49086263, 
6377420.97338432, 6375324.52698281, 6373228.15393418, 6371131.85206164, 
6369035.62364959, 6366939.46650363, 6364843.38290776, 6362747.37068454, 
6360651.43210902)), class = "data.frame", row.names = c(NA, -170L
))

Then I used the sf_multipoint {sfheaders} to convert the points to sf object :

pts_t = sf_multipoint(pts)

but the plotting does not work :

ggplot() +
     theme_void() +
     geom_sf(data = sh, fill = "grey") +
     geom_sf(data = nc_l, color = "red") +
     geom_sf(data = pts_t, color = "blue")
Error in st_transform.sfc(X[[i]], ...) : 
  cannot transform sfc object with missing crs

I'm really confused what is the best way to do this !


Solution

  • You can just use geom_point for the pts data frame:

    ggplot() +
        theme_void() +
        geom_sf(data = sh, fill = "grey") +
        geom_point(data = pts, aes(x, y), color = "red")
    

    enter image description here

    If you really want to convert pts to an sf object so you can use geom_sf, you can do:

    library(sf)
    
    pts <- st_multipoint(as.matrix(pts)) |>
      st_sfc(crs = st_crs(sh)) |>
      st_as_sf()
    
    ggplot() +
      theme_bw() +
      geom_sf(data = sh, fill = "grey") +
      geom_sf(data = nc_l, color = "red") +
      geom_sf(data = pts, color = "blue")
    

    enter image description here