用R语言从Radolan数据中制作GPS经纬度查找表，将立体投影转换为R中的经纬度网格

Question

我使用这个例子在R中加载Radolan数据。

例子, RADOLAN规范

 binary_filepath <- paste0("raa01-rw_10000-1708091950-dwd---bin")
    header_end <- regexpr("\003", readLines(binary_filepath, 1))

    rb_stream <- file(binary_filepath, "rb")
    skip_temp <- readBin(rb_stream, "raw", n = header_end, endian = "little")
    rb_data <- readBin(rb_stream, "raw", size = 1, n = 900*900*2, endian = "big")
    rb_data <- rawToBits(rb_data)

    rbi <- lapply(seq(0, 900*900-1, 1), function(x){
      sum(as.integer(rb_data[(x*16+1):(x*16+12)]) * 2**seq(0, 11, 1))
    })
    rbi <- unlist(rbi)
    rbi[1:16]
    rb <- raster(t(matrix(rbi, ncol = 900, nrow = 900)))
    rb <- flip(rb, "y")
    radolan_proj <- 
      CRS("+proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 +k=0.93301270189 +x_0=0 +y_0=0 +a=6370040 +b=6370040 +to_meter=1000 +no_defs")
    extent(rb) <- extent(-523.4622, 376.5378, -4658.645, -3758.645)
    projection(rb) <- radolan_proj

RADOLAN数据是一种立体投影网格。该网格中的每个点间隔1公里。

class       : RasterLayer 
dimensions  : 900, 900, 810000  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : -523.4622, 376.5378, -4658.645, -3758.645  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=stere +lat_0=90 +lat_ts=90 +lon_0=10 +k=0.93301270189 +x_0=0 +y_0=0 +a=6370040 +b=6370040 +to_meter=1000 +no_defs 
data source : in memory
names       : layer 
values      : 0, 2500  (min, max)

我希望为GPS数据制作一个纬度/经度查找表。如何将此rb对象转换为纬度/经度网格？

谢谢！

- Andreas

Answer 1

这个rb对象是一个栅格。如果您想要使用非投影的地理坐标（wgs84），您只需要重新投影它：

library(rgdal)
library(raster)
library(maps) # to add country limits 

# raster with Lat/Long non projected coordinates
rb_wgs84 <- projectRaster(rb, crs = "+proj=longlat +datum=WGS84 +no_defs")
# Show the raster in a plot
plot(rb_wgs84)
maps::map("world", add = TRUE)

# Get coordinates and values in a dataframe
xyz <- data.frame(coordinates(rb_wgs84), z = values(rb_wgs84))

这是您正在寻找的吗？

- Sébastien Rochette

---

是的，我正在寻找这种转换。但不知何故，数值丢失了或者我漏掉了什么？> rb_wgs84 [1:16] [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

rbi [1:16] [1] 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 / 您创建的框架中的Z值始终为NA。我想能够在数据框中搜索GPS坐标并读取该值。

- Andreas

光栅（raster）是一个矩形矩阵。如果您查看WGS84中的投影光栅，您会发现其形状并不呈矩形。这是由于立体视图投影的缘故。一些区域在矩形WGS84中没有值，因为这些是新位置。有关投影的更多信息，请参见https://en.wikipedia.org/wiki/Map_projection。因此，在最终的光栅中有NA值，但正如该图所示，还有0和2500值。您可以在数据框中选择它们：`head(xyz[!is.na(xyz$z),])`。 - Sébastien Rochette

只是一个额外的提示，但我不确定你想要实现什么，你可以使用函数raster::extract来获取特定地理位置的z值。如果你有一些来自外部数据集的GPS坐标，可以在rb_wgs84上使用它：raster::extract(rb_wgs84, mydata)。 - Sébastien Rochette

谢谢。我想要能够看到特定坐标处的雷达值。例如，我的纬度是48.77582，经度是9.182987。我希望能够看到雷达在该坐标处产生的数值。 - Andreas