代码之家  ›  专栏  ›  技术社区  ›  Andreas

从R中的Radolan数据制作GPS Lat/Long查找表,将赤平投影转换为R中的Lat长网格

  •  4
  • Andreas  · 技术社区  · 8 年前

    im使用此示例在R中加载Radolan数据

    Example , RADOLAN Specification

     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
    

    该网格中的每个点相距1公里。

     1 x 1 km-Raster

    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)
    

    如何将rb对象转换为Lat/Long栅格?

    非常感谢。

    1 回复  |  直到 7 年前
        1
  •  3
  •   Sébastien Rochette    8 年前

    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)
    

    Map in wgs84

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

    Map data

    这就是你要找的吗?