R语言 如何使用地形包为栅格数据指定正确的投影和范围?

wztqucjr  于 2022-12-06  发布在  其他
关注(0)|答案(2)|浏览(278)

我尝试使用terra r包读取tif文件,使用以下代码

hh <- rast("imagery_HH.tif")
#> Warning message:
#> [rast] unknown extent 
hh
#> class       : SpatRaster 
#> dimensions  : 8371, 8946, 1  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : 0, 8946, 0, 8371  (xmin, xmax, ymin, ymax)
#> coord. ref. :  
#> source      : imagery_HH.tif 
#> name        : imagery_HH

使用函数terra::describe("imagery_HH.tif"),我得到了以下信息:

[4] "Size is 8946, 8371"                                                      
   [5] "GCP Projection = "                                                       
   [6] "GEOGCRS[\"WGS 84\","                                                     
   [7] "    DATUM[\"World Geodetic System 1984\","                               
   [8] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"                     
   [9] "            LENGTHUNIT[\"metre\",1]]],"                                  
  [10] "    PRIMEM[\"Greenwich\",0,"                                             
  [11] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"                      
  [12] "    CS[ellipsoidal,2],"                                                  
  [13] "        AXIS[\"geodetic latitude (Lat)\",north,"                         
  [14] "            ORDER[1],"                                                   
  [15] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                  
  [16] "        AXIS[\"geodetic longitude (Lon)\",east,"                         
  [17] "            ORDER[2],"                                                   
  [18] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"                  
  [19] "    USAGE["                                                              
  [20] "        SCOPE[\"Horizontal component of 3D system.\"],"                  
  [21] "        AREA[\"World.\"],"                                               
  [22] "        BBOX[-90,-180,90,180]],"                                         
  [23] "    ID[\"EPSG\",4326]]"                                                  
  [24] "Data axis to CRS axis mapping: 2,1"

如果仔细观察,我们会发现缺少坐标参考,分辨率显示为1 x 1,范围不正确。但如果在QGIS中打开tif文件,则会显示以下crs为EPSG的属性:4326

现在如何使用terra R软件包以正确的坐标系、分辨率和范围读取tif文件?
这里是下载raster file的链接

jq6vz3qz

jq6vz3qz1#

下面的内容表明这不是一个常规的栅格。它显示了GCPs(特定栅格像元的坐标),如果需要这些,您可能没有矩形范围或恒定分辨率。(我没有检查它们是否是,但您可以)。
阅读这样的文件需要一种与读取常规光栅文件不同的方法,这是我第一次看到这样的文件,“terra”目前不支持;我会把它列在待办事项清单上。

terra::describe("imagery_HH.tif")[31:40]
 [1] "Data axis to CRS axis mapping: 2,1"                  
 [2] "GCP[  0]: Id=1, Info="                               
 [3] "          (0,0) -> (78.591314,29.400624,0)"          
 [4] "GCP[  1]: Id=2, Info="                               
 [5] "          (357.84,0) -> (78.52592634,29.41112936,0)" 
 [6] "GCP[  2]: Id=3, Info="                               
 [7] "          (715.68,0) -> (78.4607346,29.4215638,0)"   
 [8] "GCP[  3]: Id=4, Info="                               
 [9] "          (1073.52,0) -> (78.39539736,29.43198708,0)"
[10] "GCP[  4]: Id=5, Info="

(and如此直到行1383)。

6rvt4ljy

6rvt4ljy2#

正如@bretauv所指出的,投影加载正确。我不知道为什么terra::rast不能正确读取范围,但如果可以安慰一下的话,raster::raster不能读取它,所以这不仅仅是terra的问题。无论如何,如果你知道范围,可以手动设置范围,在本例中你知道。

hh <- terra::rast("imagery_HH.tif")
terra::set.ext(
  x = hh, 
  value = c(76.6811227745188262,
            78.59105666365414556,
            27.9827663027027924,
            29.6529629093873979)
)

相关问题