0

I am working in R. I need to geocode a raster image having a layer of latitude and a layer of longitude (that i’m importing with “raster” package). To be more precise, I have the first layer with an image (for example 2000x2000) that is not geocoded. I also have two more layers (raster layers) with the same dimensions (2000x2000) of the the first layer (the image) but one of these just cointains latitude and the other one just contains longitude. I would like to use latitude and longitude layers to georeference the image layer.

In few words I need to do in R what in ENVI is called GLT (Geographic Lookup Table). Can please someone help me?

um675
  • 3
  • 2
  • If you've got latitude and longitude I don't see what you have to do to "geocode" it? Can you perhaps show us some more information about your raster data? What does GLT do, precisely? – Spacedman Oct 15 '20 at 13:28
  • If I understand correctly you need to add a two bands to a raster which contain the latitude and longitude? – JonasV Oct 15 '20 at 14:52
  • Thanks for answering me. I have the first layer with an image (for example 2000x2000) that is not geocoded. I also have two more layers (raster layers) with the same dimension(2000x2000) of the the first layer (the image) but one of this just cointains latitude and the other one contains longiutude. I would like to use latitude and longitude layers to georeference the image layer. GLT works like this: -Give me the image layer -Give me latitude layer -Give me longitude layer And it gives you back a geocoded image like magic :)

    Thank you for your time

     – um675 Oct 15 '20 at 15:25
  • Is it something like the opposite of raster::xyFromCell()? – Elio Diaz Oct 15 '20 at 15:45
  • This might help. https://gis.stackexchange.com/questions/23841/create-a-raster-with-georeferenced-information-in-r. Maybe try grabbing the minimum/maximum values from both the latitude/longitude rasters and use them to set the extent – Nick Oct 15 '20 at 16:13
  • Are the latitudes and longitudes regular grids in the raster? In other words are all the values in any row of the longitude grid the same? And all the values in the columns of the latitude grid? Otherwise its possible you've got a non-rectilinear grid. What does plot(values(long), values(lat)) look like (if those are the lat-long raster)? A grid? Edit your question and put the plot there. – Spacedman Oct 15 '20 at 17:43
  • https://gis.stackexchange.com/questions/103116/map-project-a-raster-having-separate-latitude-and-longitude-raster-bands I think this will help you – Elio Diaz Oct 15 '20 at 18:51
  • The "stars" package has some support for non rectilinear rasters if that's what we've got here. Often though a non-rectilinear lat-long set of coordinates is rectilinear in some transformed coordinates, and if you can transform back then the geocoding is trivial. – Spacedman Oct 15 '20 at 19:24
  • Thank you to you all @Spacedman, so I confirm it’s regular gridded now I can’t upload image you requested me but I confirm it’s regular gridded ( lat[i,] it’s always the same value lon[,i] it’s always the vale value). – um675 Oct 16 '20 at 11:52
  • @Nick,kind of works, but it’s not that precise. It doesn’t really match google earth base map). I also need to rotate and flip the image first. – um675 Oct 16 '20 at 11:55

1 Answers1

0

First lets make some test data. A 3x4 raster with no coordinate system. R will place this on a unit square. We'll also fake up some lat-long rasters:

r = raster(matrix(1:12, 3,4))
latitude = raster(matrix(rep(c(seq(34, 37,len=3)),4),ncol=4))
longitude = raster(matrix(rep(seq(10,23,length=4),3),ncol=4, byrow=TRUE))

enter image description here

Now we get the values from longitude and latitude:

 vlong = values(longitude)
 vlat = values(latitude)

and work out the cell size. In longitude that's the difference between the first two values in the raster, and in latitude its the difference between the values in the first and the "1 + number of columns" value:

 xcell = diff(vlong[1:2])
 ycell = diff(vlat[c(1,1+ncol(latitude))])
 xcell
# [1] 4.333333
 ycell
# [1] 1.5

Now set the extent of r to the min/max lat/long values plus or minus half a cell size in the right direction:

 extent(r) = c(min(vlong)-xcell/2, max(vlong)+xcell/2, min(vlat)-ycell/2, max(vlat)+ycell/2)
 print(r)
# class      : RasterLayer 
# dimensions : 3, 4, 12  (nrow, ncol, ncell)
# resolution : 4.333333, 1.5  (x, y)
# extent     : 7.833333, 25.16667, 33.25, 37.75  (xmin, xmax, ymin, ymax)

Now we can test this by plotting the raster then adding the points referred to by corresponding points in the lat-long rasters:

> plot(r)
> points(longitude[], latitude[])

enter image description here

which looks dead on. If this doesn't work for you then possibly your lat and long values aren't equally spaced.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • Thank you very much! But what if I write crs(r) or if I want to export it as geotiff writeRaster(x,”nameraster.tif”,format=“GTiff”) ? Thank you very much again! – um675 Oct 18 '20 at 12:34
  • Yes, if you set the CRS with crs(r) = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" then you should be able to write it as a GeoTIFF. – Spacedman Oct 18 '20 at 14:39
  • Thank you very much! Now I don’t have my computer with me. Next week I will try this solution, but it looks like the complete solution! As soon as possible I will let you know! – um675 Oct 18 '20 at 15:14
  • Both your solution and the solution proposed by @Nick work. There is still a problem regarding the fact that I need to flipLR and rotate -90 the image before geocoding, but if the footprint is not a square footprint it is impossible to rotate the image and not rotating the lat and lon raster. I have no idea why I need to flip and rotate the image, any clue? – um675 Oct 27 '20 at 18:31