多边形的非唯一属性的中位栅格值

huangapple go评论60阅读模式
英文:

median raster value for polygons by attribute that is not unique

问题

给定这个栅格和多边形数据集:

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
fs <- system.file("ex/lux.shp", package="terra")
v <- vect(fs)

如何提取多边形 v 中具有相同 "NAME_1" 值的栅格单元的中位数值,并将其添加到 v 的属性中?

请注意,"NAME_1" 在 v 中不是唯一的:

v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 12, 6  (geometries, attributes)
# extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
# source      : lux.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
# type        : <num>    <chr> <num>    <chr> <num> <int>
# values      :     1 Diekirch     1 Clervaux   312 18081
#                   1 Diekirch     2 Diekirch   218 32543
#                   1 Diekirch     3  Redange   259 18664
英文:

Given this this raster and polygon dataset

 f &lt;- system.file(&quot;ex/elev.tif&quot;, package=&quot;terra&quot;) 
 r &lt;- rast(f)
 fs &lt;- system.file(&quot;ex/lux.shp&quot;, package=&quot;terra&quot;)
 v &lt;- vect(fs)

How can I extract the median value of the raster cells for polygons with the same value of "NAME_1" in v, and add this to the attributes of v?

Please note that "NAME_1" is not unique in v:

v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 12, 6  (geometries, attributes)
# extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
# source      : lux.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP
# type        : &lt;num&gt;    &lt;chr&gt; &lt;num&gt;    &lt;chr&gt; &lt;num&gt; &lt;int&gt;
# values      :     1 Diekirch     1 Clervaux   312 18081
#                   1 Diekirch     2 Diekirch   218 32543
#                   1 Diekirch     3  Redange   259 18664

答案1

得分: 1

为了计算所有几何体的中位数,您可以执行以下操作:

library(terra)
r <- rast(system.file("ex/elev.tif", package="terra"))
v <- vect(system.file("ex/lux.shp", package="terra"))

v$median = extract(r, v, fun=median, na.rm=TRUE, ID=FALSE)
v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 12, 7  (几何体, 属性)
# extent      : 5.74414, 6.528252, 49.44781, 50.18162  (最小x, 最大x, 最小y, 最大y)
# source      : lux.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP median
# type        : <num>    <chr> <num>    <chr> <num> <int>  <num>
# values      :     1 Diekirch     1 Clervaux   312 18081  467.1
#                   1 Diekirch     2 Diekirch   218 32543  333.9
#                   1 Diekirch     3  Redange   259 18664  377.4

由于"NAME_1"对于多个几何体具有相同的值,我们需要多做一些工作。您可以选择执行以下操作之一:

a <- aggregate(v, "NAME_1")
e <- extract(r, a, fun=median, na.rm=TRUE)
e$NAME_1 <- a$NAME_1[e$ID]
e$ID <- NULL
ve <- merge(v, e, by="NAME_1")

或者

ee <- extract(r, v, na.rm=TRUE)
ee$NAME_1 <- v$NAME_1[ee$ID]
ee <- aggregate(ee[,2,drop=FALSE], ee[,"NAME_1",drop=FALSE], median, na.rm=TRUE)
vee <- merge(v, ee, by="NAME_1")

values(vee)
#         NAME_1 ID_1 ID_2           NAME_2 AREA    POP elevation
#1      Diekirch    1    1         Clervaux  312  18081     422.0
#2      Diekirch    1    2         Diekirch  218  32543     422.0
#3      Diekirch    1    3          Redange  259  18664     422.0
#4      Diekirch    1    4          Vianden   76   5163     422.0
#5      Diekirch    1    5            Wiltz  263  16735     422.0
#6  Grevenmacher    2    6       Echternach  188  18899     288.5
#7  Grevenmacher    2    7           Remich  129  22366     288.5
#8  Grevenmacher    2   12     Grevenmacher  210  29828     288.5
#9    Luxembourg    3    8         Capellen  185  48187     314.0
#10   Luxembourg    3    9 Esch-sur-Alzette  251 176820     314.0
#11   Luxembourg    3   10       Luxembourg  237 182607     314.0
#12   Luxembourg    3   11           Mersch  233  32112     314.0
英文:

To compute the median for all geometries, you can do

library(terra)
r &lt;- rast(system.file(&quot;ex/elev.tif&quot;, package=&quot;terra&quot;) )
v &lt;- vect(system.file(&quot;ex/lux.shp&quot;, package=&quot;terra&quot;))

v$median = extract(r, v, fun=median, na.rm=TRUE, ID=FALSE)
v
# class       : SpatVector 
# geometry    : polygons 
# dimensions  : 12, 7  (geometries, attributes)
# extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
# source      : lux.shp
# coord. ref. : lon/lat WGS 84 (EPSG:4326) 
# names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP median
# type        : &lt;num&gt;    &lt;chr&gt; &lt;num&gt;    &lt;chr&gt; &lt;num&gt; &lt;int&gt;  &lt;num&gt;
# values      :     1 Diekirch     1 Clervaux   312 18081  467.1
#                   1 Diekirch     2 Diekirch   218 32543  333.9
#                   1 Diekirch     3  Redange   259 18664  377.4

Since "NAME_1" has the same value for multiple geometries, we have to do a bit more work. You can do either

a &lt;- aggregate(v, &quot;NAME_1&quot;)
e &lt;- extract(r, a, fun=median, na.rm=TRUE)
e$NAME_1 &lt;- a$NAME_1[e$ID]
e$ID &lt;- NULL
ve &lt;- merge(v, e, by=&quot;NAME_1&quot;)

or

ee &lt;- extract(r, v, na.rm=TRUE)
ee$NAME_1 &lt;- v$NAME_1[ee$ID]
ee &lt;- aggregate(ee[,2,drop=FALSE], ee[,&quot;NAME_1&quot;,drop=FALSE], median, na.rm=TRUE)
vee &lt;- merge(v, ee, by=&quot;NAME_1&quot;)

values(vee)
#         NAME_1 ID_1 ID_2           NAME_2 AREA    POP elevation
#1      Diekirch    1    1         Clervaux  312  18081     422.0
#2      Diekirch    1    2         Diekirch  218  32543     422.0
#3      Diekirch    1    3          Redange  259  18664     422.0
#4      Diekirch    1    4          Vianden   76   5163     422.0
#5      Diekirch    1    5            Wiltz  263  16735     422.0
#6  Grevenmacher    2    6       Echternach  188  18899     288.5
#7  Grevenmacher    2    7           Remich  129  22366     288.5
#8  Grevenmacher    2   12     Grevenmacher  210  29828     288.5
#9    Luxembourg    3    8         Capellen  185  48187     314.0
#10   Luxembourg    3    9 Esch-sur-Alzette  251 176820     314.0
#11   Luxembourg    3   10       Luxembourg  237 182607     314.0
#12   Luxembourg    3   11           Mersch  233  32112     314.0

huangapple
  • 本文由 发表于 2023年4月11日 01:52:36
  • 转载请务必保留本文链接:https://go.coder-hub.com/75979459.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定