how to properly project a rast object if extent information is missing?

Question

I'm trying to read a HDF5 file with terra, but the extent of the grid could not be read.

> rst <- terra::rast("temp.h5")
    Warning message:
    [rast] unknown extent
     
> rst
    class       : SpatRaster 
    dimensions  : 765, 700, 2  (nrow, ncol, nlyr)
    resolution  : 0.001428571, 0.00130719  (x, y)
    extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
    coord. ref. : lon/lat WGS 84 
    sources     : temp.h5://image_data  
                  temp.h5://image_data  
    varnames    : image_data 
                  image_data 
    names       : image_data, image_data

with the function terra::describe("temp.h5") I obtained the information below:

 [4] "Metadata:"                                                                                                                       
 [5] "  geographic_geo_column_offset=0 "                                                                                               
 [6] "  geographic_geo_dim_pixel=KM,KM"                                                                                                
 [7] "  geographic_geo_number_columns=700 "                                                                                            
 [8] "  geographic_geo_number_rows=765 "                                                                                               
 [9] "  geographic_geo_par_pixel=X,Y"                                                                                                  
[10] "  geographic_geo_pixel_def=LU"                                                                                                   
[11] "  geographic_geo_pixel_size_x=1.0000013 "                                                                                        
[12] "  geographic_geo_pixel_size_y=-1.0000055 "                                                                                       
[13] "  geographic_geo_product_corners=0 49.362053 0 55.973602 10.856429 55.388977 9.0092793 48.895298 "                               
[14] "  geographic_geo_row_offset=3649.9792 "                                                                                          
[15] "  geographic_map_projection_projection_indication=Y"                                                                             
[16] "  geographic_map_projection_projection_name=STEREOGRAPHIC"                                                                       
[17] "  geographic_map_projection_projection_proj4_params=+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0"

From the he corners of the map (line 13), I read that xmin=0, xmax=10.856429, ymin=49.362053 and ymax=55.973602. However, when setting the extent of the raster and plotting, I see that the corners of the raster do not correspond with the corners as listed in line 13 of the output above.

ext(rst) <- c(0,10.856429,49.362053,55.973602)
plot(rst)

Clearly something is going wrong here. How do I fix this?

Thanks in advance!

score 1 · Answer 1 · answered Apr 02 '22 at 18:16

The metadata show that the CRS is stereographic, this not lon/lat. We can set it like this (note that I added the "+units=km" based on the values of a and b

p <- "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.14 +b=6356.75 +x_0=0 y_0=0 +units=km"
crs(rst) <- p

Now you need to find the extent in that projection. Perhaps that is available in the file, but we can also try to guess it from the lon/lat bounding box.

m <- matrix(c(0, 49.362053, 0, 55.973602, 10.856429, 55.388977, 9.0092793, 48.895298), ncol=2, byrow=T)    
v <- vect(m, crs="+proj=longlat +datum=WGS84")
# I had to change the a and b parameters to m to make this work
p <- "+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378137 +b=6356750 +x_0=0 y_0=0"
s <- project(v, p)
e <- extent(s)
e 
#SpatExtent : 0, 700000.402146425, -4415000.95326333, -3649996.04313886 (xmin, xmax, ymin, ymax)

And now you could try

ext(rst) <- as.vector(e)/1000

And you may also need to consider this:

geographic_geo_row_offset=3649.9792 "

Overall, it would seem best to go to the source and ask there.

how to properly project a rast object if extent information is missing?

1 Answers1