我有一个栅格，我想将其值更改为每公顷像素值，其中面积来自像素的大小。因此，如果一个像素是 10 公顷，其当前值为 100 个任意单位，那么现在它将是 10 个任意单位/公顷。
栅格示例：
r <- rast(nrows=18, ncols=36)
v <- 1:ncell(r)
v[200:400] <- NA
values(r) <- v
此函数应采用栅格并乘以 1/ha。
px2ha <- function(raster) {
raster * 10000/prod(res(raster))
}
那么像这样使用 cellSize 怎么样？
r2/cellSize(r2, unit = 'ha')
我没有得到相同的答案。
r2 = project(x = r, y = 'epsg:2950') # Put raster in meters
plot((r2))
plot(px2ha(r2))
plot(r2/cellSize(r2, unit = 'ha'))
可以使用 cellSize 获取每面积（公顷）的像素值吗？
我有此代码尝试从此处找到的法律中提取物种https://laws.justice.gc.ca/fra/lois/S-15.3/TexteComplet.html
但是，我无法让 html_nodes 找到每个部分
section <- div_content %>% html_nodes(xpath = paste0("//h2[contains(text(), '", header, "')]/following-sibling::div[contains(@class, 'ProvisionList')]"))
基本上，我找不到获取文本内容并匹配其他部分的方法。我尝试添加“
”标签并查找每个部分的文本，但没有成功（获取{xml_nodeset (0)}）
我尝试获取 id 为“425426”的 div 中的数据，然后在 scheduleLabel 中从 scheduleTitleText 中获取文本。我需要另一列用于 SchedHeadL1（即包含物种的部分的标题）和 BilingualGroupTitleText 中的文本（说明动物或植物的组...）。然后提供一个嵌套的物种列表（这里我将物种与法语名称、拉丁语和英语分开）
library(rvest)
library(dplyr)
library(stringr)
# URL of the webpage
url <- "https://laws.justice.gc.ca/fra/lois/S-15.3/TexteComplet.html"
# Read the webpage content
webpage <- read_html(url)
# Extract the div with id "425426"
div_content <- webpage %>% html_node("#425426")
# Extract the header h2 with class "scheduleTitleText" from the class "scheduleLabel" and id "h-425427"
schedule_label <- div_content %>% html_node("h2.scheduleLabel#h-425427") %>% html_text()
# Extract all h2 headers with class "SchedHeadL1"
headers <- div_content %>% html_nodes("h2.SchedHeadL1") %>% html_text()
# Use str_extract to extract the "PARTIE #" part
partie_numbers <- str_extract(headers, "PARTIE \\d+")
# Use str_remove to remove the "PARTIE #" part from the original strings
descriptions <- str_remove(headers, "PARTIE \\d+")
# Combine into a data frame
result <- data.frame(Partie = partie_numbers, Description = descriptions, stringsAsFactors = FALSE)
headers_prep = result |>
unite(pd, Partie, Description, sep = "<br>") |> pull(pd)
# Initialize lists to store the extracted data
group_titles <- list()
item_first <- list()
item_second <- list()
scientific_names <- list()
latin_names <- list()
# Loop through each header to extract the associated content
for (header in headers) {
# Extract the section associated with the current header
section <- div_content %>% html_nodes(xpath = paste0("//h2[contains(text(), '", header, "')]/following-sibling::div[contains(@class, 'ProvisionList')]"))
# Extract BilingualGroupTitleText within the section
group_title <- section %>% html_nodes(".BilingualGroupTitleText") %>% html_text()
group_titles <- c(group_titles, group_title)
# Extract BilingualItemFirst within the section
item_first_section <- section %>% html_nodes(".BilingualItemFirst") %>% html_text()
item_first <- c(item_first, item_first_section)
# Extract BilingualItemSecond within the section
item_second_section <- section %>% html_nodes(".BilingualItemSecond") %>% html_text()
item_second <- c(item_second, item_second_section)
# Extract otherLang (scientific names) within the section
scientific_name_section <- section %>% html_nodes(".otherLang") %>% html_text()
scientific_names <- c(scientific_names, scientific_name_section)
# Extract scientific Latin names from BilingualItemFirst
latin_name_section <- str_extract(item_first_section, "\\(([^)]+)\\)") %>% str_replace_all("[()]", "")
latin_names <- c(latin_names, latin_name_section)
}
# Ensure all columns have the same length by repeating the last element if necessary
max_length <- max(length(headers), length(group_titles), length(item_first), length(item_second), length(scientific_names), length(latin_names))
schedule_label <- rep(schedule_label, length.out = max_length)
headers <- rep(headers, length.out = max_length)
group_titles <- rep(group_titles, length.out = max_length)
item_first <- rep(item_first, length.out = max_length)
item_second <- rep(item_second, length.out = max_length)
scientific_names <- rep(scientific_names, length.out = max_length)
latin_names <- rep(latin_names, length.out = max_length)
# Create a data frame
data <- data.frame(
ScheduleLabel = schedule_label,
Header = headers,
GroupTitle = group_titles,
ItemFirst = item_first,
ItemSecond = item_second,
ScientificName = scientific_names,
LatinName = latin_names,
stringsAsFactors = FALSE
)
我有 2 个列表：
a = list(data.frame(a = 1:10, b = 21:30), data.frame(a = 11:20, b = 31:40))
b = list(data.frame(c = letters[1:10]), data.frame(c = LETTERS[1:10]))
我想自动对它们进行 cbind，以便 c 列与 a 列表中的每个元素配对。
这是我的期望（手动）
c = list(data.frame(a = 1:10, b = 21:30, c = letters[1:10]), data.frame(a = 11:20, b = 31:40, c = LETTERS[1:10]))
当尝试这个时，
do.call(cbind, a, b)
我明白了
Error in if (quote) args <- lapply(args, enquote) :
the condition has length > 1
预期输出：
c
[[1]]
a b c
1 1 21 a
2 2 22 b
3 3 23 c
4 4 24 d
5 5 25 e
6 6 26 f
7 7 27 g
8 8 28 h
9 9 29 i
10 10 30 j
[[2]]
a b c
1 11 31 A
2 12 32 B
3 13 33 C
4 14 34 D
5 15 35 E
6 16 36 F
7 17 37 G
8 18 38 H
9 19 39 I
10 20 40 J
将 sf 对象打印到控制台时，有一个带有Geodetic CRS或 的标题Porjected CRS，其中包含 CRS 的良好格式名称。有没有办法提取此信息，而不是
p1 = st_point(c(7,52))
p2 = st_point(c(-30,20))
sfc = st_sfc(p1, p2, crs = 4326)
sfc
Geometry set for 2 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -30 ymin: 20 xmax: 7 ymax: 52
Geodetic CRS: WGS 84 <- this line
POINT (7 52)
POINT (-30 20)
st_transform(sfc, 3857)
Geometry set for 2 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -3339585 ymin: 2273031 xmax: 779236.4 ymax: 6800125
Projected CRS: WGS 84 / Pseudo-Mercator <- this line
POINT (779236.4 6800125)
POINT (-3339585 2273031)
nc = st_read(system.file("shape/nc.shp", package="sf"))
nc
Simple feature collection with 100 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS: NAD27 <- this line
First 10 features:
[...]
看看如何st_crs(sfc)$input给予"EPSG:4326"而不是WGS 84。然而，st_crs(nc)$input给予"NAD27"，但st_crs(st_transform(sfc, 3857))$input给予 "EPSG:3857"而不是Projected CRS: WGS 84 / Pseudo-Mercator
地球上是否也存在与栅格数据相当的东西？
我有一个代码，其中有一个场景（sc列）和一个代码名称（cd列）。我想计算所有场景与第一个场景之间的差异（所以场景2-1，3-1等等。我还有多个场景......）。有没有一种快速的方法来编码将要比较哪个场景，例如，如果我现在所有其他场景的比较场景现在是场景 2，那将是1-2，3-2等等？
但我也想确保在场景之间只减去相同的代号名称（cd），所以（sc1.cd1 - sc2.cd1），（sc1.cd2 - sc2.cd2）等等。
sc我还想要了解每个代码中 1 与其他场景之间的百分比差异。
目前，我的代码假设当我比较场景时，它们cd的顺序相同。有没有办法将场景和代码一起比较，而不是使用数据子集？
set.seed(123456)
vals = as.vector(mapply(FUN = rpois, n = 10, lambda = c(5, 10, 100)))
df.test = data.frame(cd = rep(1:10,3), sc = rep(1:3, each =10 ), vals)
# New empty column
df.test$delta = NA
# Subset each scenarios and subtract for the first scenario
df.test[df.test$sc ==1,'delta'] = df.test[df.test$sc ==1,'vals']-df.test[df.test$sc ==1,'vals']
df.test[df.test$sc ==2,'delta'] = df.test[df.test$sc ==2,'vals']-df.test[df.test$sc ==1,'vals']
df.test[df.test$sc ==3,'delta'] = df.test[df.test$sc ==3,'vals']-df.test[df.test$sc ==1,'vals']
# Get difference
df.test$diff = abs(abs(df.test$delta) - df.test$vals)
# Calculate percentage
df.test$delta.perc = abs(df.test$delta) / df.test$diff +1
df.test$diff * df.test$delta.perc
# Percentage increase
df.test$delta.perc.increase = df.test$delta.perc *100
我正在使用 mapview 包构建地图。但是，当我使用 st_sample 并生成点时，有些点没有被映射。另外，st_buffer合并了一些区域（它们应该是分开的）。
library(sf)
library(mapview)
library(magrittr)
set.seed(145)
chemins.32188 = st_zm(chemins2) %>% st_transform(crs = 32188)
rd.pts.sentiers <- st_sample(x = chemins.32188, size = 10)
buff.pt = st_buffer(x = rd.pts.sentiers,
dist = units::set_units(x = 100,
value = "m"))
mapviewOptions(fgb = FALSE)
mapview(chemins.32188) +
mapview(rd.pts.sentiers) +
mapview(buff.pt)
当我使用绘图时，它可以工作......但缓冲区仍然被合并。
plot(chemins.32188$geom)
plot(buff.pt, add= TRUE, col = scales::alpha("blue", .5))
plot(rd.pts.sentiers, add= TRUE, pch = 19, col = "red")
这是数据
chemins2 = structure(list(Name = "Path Measure", description = NA_character_,
timestamp = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), begin = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), end = structure(NA_real_, class = c("POSIXct", "POSIXt"
)), altitudeMode = NA_character_, tessellate = -1L, extrude = 0L,
visibility = -1L, drawOrder = NA_integer_, icon = NA_character_,
sentier = NA_character_, layer = "path Parc — Path Measure",
path = "~/Desktop/path Parc.kml|layername=Path Measure|geometrytype=LineString25D|uniqueGeometryType=yes",
geom = structure(list(structure(list(structure(c(299434.130314458,
300531.43655987, 299434.130314458, 5046816.92873926, 5047257.53966348,
5046816.92873926), dim = 3:2), structure(c(299434.130314458,
300401.157828111, 299434.130314458, 5046816.92873926, 5046249.12770288,
5046816.92873926), dim = 3:2)), class = c("XY", "MULTILINESTRING",
"sfg"))), n_empty = 0L, crs = structure(list(input = "EPSG:32188",
wkt = "PROJCRS[\"NAD83 / MTM zone 8\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"MTM zone 8\",\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\",-73.5,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9999,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",304800,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"easting (E(X))\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"northing (N(Y))\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Canada - Quebec between 75°W and 72°W.; Canada - Ontario - east of 75°W.\"],\n BBOX[44.98,-75,62.53,-72]],\n ID[\"EPSG\",32188]]"), class = "crs"), class = c("sfc_MULTILINESTRING",
"sfc"), precision = 0, bbox = structure(c(xmin = 299434.130314458,
ymin = 5046249.12770288, xmax = 300531.43655987, ymax = 5047257.53966348
), class = "bbox"))), row.names = 1L, sf_column = "geom", agr = structure(c(Name = NA_integer_,
description = NA_integer_, timestamp = NA_integer_, begin = NA_integer_,
end = NA_integer_, altitudeMode = NA_integer_, tessellate = NA_integer_,
extrude = NA_integer_, visibility = NA_integer_, drawOrder = NA_integer_,
icon = NA_integer_, sentier = NA_integer_, layer = NA_integer_,
path = NA_integer_), class = "factor", levels = c("constant",
"aggregate", "identity")), class = c("sf", "data.frame"))
我有一个包含命名数据框的命名数据列表。
mlst = list(test = list(df = data.frame(x = 1:10, y = 21:30), list(c(1:2))),
test2 = list(df = data.frame(x = 1:10, y = 21:30), list(c(1:2))))
mlst$test$df
mlst$test2$df
该功能do.call(rbind,mlst)不起作用，因为它不仅仅拾取数据帧。如何在没有 for 循环的情况下仅绑定数据帧的行？
在这个网站上https://www.quebec.ca/agriculture-environnement-et-ressources-naturelles/faune/gestion-faune-habitats-fauniques/especes-fauniques-menacees-vulnerables/liste，有我的物种表想要提取。
library(rvest)
sp.list = "https://www.quebec.ca/agriculture-environnement-et-ressources-naturelles/faune/gestion-faune-habitats-fauniques/especes-fauniques-menacees-vulnerables/liste"
# Get website
wp.list = read_html(species.list)
# Extract name of sections
headers = wp.list %>% html_elements("h3") %>% html_text2() %>% .[1:24]
# Get tables
tab = read_html(sp.list) %>% html_table(header = TRUE)
# Name tables
names(tab) = headers
# Combine tables
tab.gr = dplyr::bind_rows(tab, .id = "group")
这使：
tab.gr
# A tibble: 180 × 3
group Espèce `Nom latin`
<chr> <chr> <chr>
1 Mollusques Anodonte du gaspareau Utterbackiana implicata
2 Mollusques Obovarie olivâtre Obovaria olivaria
3 Insectes Bourdon à tache rousse Bombus affinis
4 Insectes Coccinelle à neuf points Coccinella novemnotata
我能够获取节标题h2，但无法将它们与每个h3节相关联
get.section = wp.list %>% html_nodes('.frame, .frame-default, .frame-type-textmedia, .frame-layout-0')
pas.dans.cette.page = !grepl(pattern = "Dans cette", x = get.section)
subset.listes = get.section[pas.dans.cette.page]
sections.tables = subset.listes[grep(pattern = "Liste des esp", x = subset.listes)]
sections.tables %>% html_elements("h2") %>% html_text2()
[1] "Liste des espèces menacées"
[2] "Liste des espèces vulnérables"
[3] "Liste des espèces susceptibles d’être désignées comme menacées ou vulnérables"
那么我怎样才能获得标题（例如，“Liste des espèces menacées”）及其组（例如，“Mollusques”）及其表格？
我有列表（元素数量不均匀）：
my.list = list(col1 = c("CC", "CT", "TT"),
col2 = c("GG", "GT"),
col3 = c("CC", "CT"),
col4 = c("CC", "CG", "GG"),
col5 = c("AC", "CC"),
col6 = "GG")
$col1
[1] "CC" "CT" "TT"
$col2
[1] "GG" "GT"
$col3
[1] "CC" "CT"
$col4
[1] "CC" "CG" "GG"
$col5
[1] "AC" "CC"
$col6
[1] "GG"
可以将其转换为数据框：
mylist.df = plyr::ldply(my.list, rbind)
names(mylist.df) <- c("cols","g1", "g2", "g3")
我想使用mylistor对下面的数据框进行子集化mylist.df。基本上，我想保留每个值至少有一个元素的每一行mylist：
df.to.subset = structure(list(IDs = c("ID1", "ID2", "ID3", "ID4", "ID5", "ID6"),
gr = c("gr1", "gr1", "gr1", "gr1", "gr1", "gr1"),
var = c(-3.451, -3.469, -3.837, -3.344, -3.904, -3.943),
col1 = structure(c(1L, 2L, 3L, 1L, 2L, 2L), levels = c("CC", "CT", "TT"), class = "factor"),
col2 = structure(c(1L, 1L, 2L, 3L, 3L, 3L), levels = c("GG", "GT", "TT"), class = "factor"),
col3 = structure(c(1L, 2L, 1L, 1L, 1L, 1L), levels = c("CC", "CT"), class = "factor"),
col4 = structure(c(2L, 2L, 2L, 2L, 2L, 2L), levels = c("CC", "CG", "GG"), class = "factor"),
col5 = structure(c(1L, 2L, 2L, 2L, 2L, 2L), levels = c("AC", "CC"), class = "factor"),
col6 = structure(c(1L, 1L, 2L, 1L, 1L, 1L), levels = c("GG","AA"), class = "factor")),
row.names = c(NA,
-6L), class = c("tbl_df", "tbl", "data.frame"))
IDs gr var col1 col2 col3 col4 col5 col6
<chr> <chr> <dbl> <fct> <fct> <fct> <fct> <fct> <fct>
1 ID1 gr1 -3.45 CC GG CC CG AC GG
2 ID2 gr1 -3.47 CT GG CT CG CC GG
3 ID3 gr1 -3.84 TT GT CC CG CC AA
4 ID4 gr1 -3.34 CC TT CC CG CC GG
5 ID5 gr1 -3.90 CT TT CC CG CC GG
6 ID6 gr1 -3.94 CT TT CC CG CC GG
（最终结果将是）
IDs gr var col1 col2 col3 col4 col5 col6
ID1 gr1 -3.45 CC GG CC CG AC GG
ID2 gr1 -3.47 CT GG CT CG CC GG
另外，我想重新调整每一列的级别以df.to.subset匹配此数据框中的级别：
factor.levels.cols = structure(list(cols = c("col1", "col2", "col3", "col4", "col5", "col6"),
g1 = c("CC", "GG", "CC", "CC", "AA", "AA"),
g2 = c("CT", "GT", "CT", "CG", "AC", "AG"),
g3 = c("TT", "TT", "TT", "GG", "CC", "GG")),
row.names = c(NA, 6L), class = "data.frame")
cols g1 g2 g3
1 col1 CC CT TT
2 col2 GG GT TT
3 col3 CC CT TT
4 col4 CC CG GG
5 col5 AA AC CC
6 col6 AA AG GG
这里是否强制使用 for 循环，或者有没有办法让它更快？我有超过 1,000,000 个条目需要修改。
我想使用facet_wrap 在x 轴上为不同的组绘制不同的值。但是，使用 时scale = "free_x"，它会显示每个组内的唯一值。请参阅下面的示例：
library(ggplot2)
set.seed(1234)
nb = rep(LETTERS[1:3],4)
expand.grid(LETTERS[1:3], 1:3)
df = data.frame(ps = paste(rep(c("ps1","gc1"),6)),
cr = rep(LETTERS[1:3],4),
id = rep(x = month.abb[1:4], each = 3));df
df = df[order(df$ps),]
df[df$ps == "ps1","cr"] <- sample(LETTERS[4:6], size = 6, replace = TRUE)
df.all = rbind(df,df,df,df,df,df)
df.all$val = rnorm(n = nrow(df.all))
ggplot(df.all,
aes(x = cr,
y = val)) +
geom_violin() +
facet_wrap(vars(ps, id),
scales = "free_x",
ncol = 4) +
geom_point(position = position_jitter(width = 0.1, height = 0.1),
alpha= .3)
在此处的图像中（这是上面代码的输出），我添加了缺少的级别。基本上，有一个组“gc1”，我希望其中的级别不要下降。但我也不希望“ps1”中的级别出现在“gc1”中。所以最后 gc1 应该只显示级别 A、B、C，而 ps1 应该只显示 D、E 和 F。有没有办法做到这一点？
