CODEWITHSUNDEEP
rterra
UseR10085
UseR10085

Reputation: 8166

How to assign correct projection and extent to a raster using terra r package?

I am trying to read a tif file using terra r package using the following code

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

Using the function terra::describe("imagery_HH.tif"), I got the following information:

   [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"

If we look closely, we can see that the coordinate reference is missing and the resolution is showing 1 x 1 with the incorrect extent. But if we open the tif file in QGIS, it shows the following properties having a crs of EPSG:4326

enter image description here

Now how to read the tif file with proper coordiante system, resolution and extent using terra R package.

Upvotes: 0

Views: 1209

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47146

The below suggests that this is not a regular raster. It shows GCPs (coordinates for particular raster cells) and if these are needed you probably do not have a rectangular extent or constant resolution. (I have not checked if they are, but you could).

Reading such a file requires a different approach than reading a regular raster file. This is the first time I see a file like this, and "terra" currently does not support it; I will put it on the to-do list.

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 so on until line 1383).

Upvotes: 1

Anthony Martinez
Anthony Martinez

Reputation: 160

As noted by @bretauv, the projection is loading correctly. I'm not sure why the extent isn't being read correctly by terra::rast, but if it's any consolation, it isn't read by raster::raster, so it's not just a terra problem. In any case, the extent can be set manually if you know the extent, which in this case you do.

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

Upvotes: 0

Related Questions