这可能是一种非常低效的方法，但这是一个开始：

library(sf)

# increment in metres.
# Smaller numbers will give you a more accurate map, but will take longer to calculate.
increment <- 10*1000

# load shapefile and convert to an equal-area projection so we can work in metres
gerShp <- st_read('gadm36_DEU_shp/gadm36_DEU_0.shp') # change this to the correct path
gerShp <- st_transform(gerShp, 3035)

# calculate total area and our 30% value
totalArea <- st_area(gerShp)
thirtyPC <- totalArea * 0.3

# Plot it
plot(gerShp, col = 'lightgrey', border = 'darkgrey', max.plot=1, reset=F)

# Find the bounding box of the feature
bbox <- st_bbox(gerShp)

thisArea <- totalArea - totalArea # zero with correct units
i <- 1

# While our subarea is less than 30%...
while (thisArea < thirtyPC) {

    # Starting at bottom, create a bounding box that is smaller than full bounding box
    thisBBox <- bbox
    thisBBox['ymax'] <- thisBBox$ymin + (increment * i)
    
    # Clip shp to this bounding box
    thisSubarea <- st_crop(gerShp, y=thisBBox)
    thisArea <- st_area(thisSubarea)
    
    print(thisArea)
    
    i <- i + 1
    
}

plot(thisSubarea, max.plot=1, add=T, col='red', border=NA)
actualPercentage <- thisArea / totalArea

在本例中，actualPercentage = 0.3011579，非常接近我们的目标值30%。

library(EBImage)
library(rgdal)
ger <- readOGR("gadm36_DEU_shp/gadm36_DEU_0.shp")

# https://stackoverflow.com/questions/16496210/rotate-a-matrix-in-r#comment-23680230
fillup.rotate = function(mat, rotations = 1) {
  for(i in seq(1:rotations)) {
    mat <- t(mat[nrow(mat):1,,drop=FALSE])
  }
  return(mat)
}

# shape file
# p 0-1 percentage
# bgcolor background fill color
# fillcolor percentage fill color
# rotations orientation
fillup <- function(shape, p = 0.5, bgcolor = "#FF0000", fillcolor = "#999999", bordercolor = "#000000", rotations = 3, width = 1000, height = 1000) {

  png("shape.png", width = width, height = height)
  par(mar=c(0,0,0,0))
  plot(shape, col=bgcolor , bg = "transparent", border = bordercolor)
  dev.off()

  image <- readImage("shape.png")
  shape.raster <- as.raster(image)

  # rotations
  # 1 top down
  # 2 left to right
  # 3 bottom up
  # 4 rigth to left

  # rotate 
  shape.raster <- fillup.rotate(shape.raster, rotations)

  # find background color 
  idx <- which(shape.raster == bgcolor)
  idx.rev <- rev(idx)

  # calc percentage
  pixel.summe <- length(idx)
  pixel.p <- pixel.summe * p

  idx.p <- idx.rev[seq(from = 1, to = pixel.p)]
  shape.raster[idx.p] <- fillcolor

  rest <- 4 - (rotations %% 4)

  # rotate back
  shape.raster <- fillup.rotate(shape.raster, rest)

  return(shape.raster)
}

ger.w <- 1000
ger.h <- 1000
colors <- c("#AAAAAA", "#333333")
ger.raster <- fillup(ger
                     , p = 0.3
                     , bgcolor = "#AAAAAA"
                     , fillcolor = "#333333"
                     , bordercolor = "#000000"
                     , rotations = 3
                     , width= ger.w
                     , height = ger.h)

plot(ger.raster)