1

I am trying to import a geotiff of 8-band satellite imagery into R using raster.

g <- raster('16SEP22191919-M1BS-057310833010_01_P001.tif' )

When I call the object it displays the extent as the the image size in pixels despite having the correct CRS info (I think):

g
class : RasterLayer
band : 1 (of 8 bands)
dimensions : 4336, 10651, 46182736 (nrow, ncol, ncell)
resolution : 1, 1 (x, y)
extent : 0, 10651, 0, 4336 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : D:\Digital Globe Imagery\nasongeotiff\057310833010_01_003\057310833010_01\057310833010_01_P001_MUL\16SEP22191919-M1BS-057310833010_01_P001.tif
names : X16SEP22191919.M1BS.057310833010_01_P001
values : 0, 65535 (min, max)

If I try to display the image I get this error:

plotRGB(g, r = 5, g = 3, b = 2, axes = TRUE, stretch = "lin")
Warning message: In couldBeLonLat(x) : raster has a longitude/latitude CRS, but coordinates do not match that

This file correctly imports into in ArcMap with the appropriate coordinates, and has this info:

  • GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433],AUTHORITY["EPSG",4326]]

The correct lat and long extents should be (-120.893178348, -120.703696534, 47.7366916019, 47.797656863), but i have had no luck defining the extent in R. Unfortunately the file is way too large (>1gb) for me to post.

Edit: There is an associated .imd file. Here are its contents: version = "28.3";
generationTime = 2017-12-06T20:55:20.000000Z;
productOrderId = "057310833010_01_P001";
productCatalogId = "A010010322340E00";
childCatalogId = "2020010322340F00";
imageDescriptor = "Basic1B";
bandId = "Multi";
panSharpenAlgorithm = "None";
numRows = 4336; numColumns = 10651; productLevel = "LV1B";
productType = "Basic";
numberOfLooks = 1;
radiometricLevel = "Corrected"; bitsPerPixel = 16;
compressionType = "None";
BEGIN_GROUP = BAND_C
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 1.432870e-02; effectiveBandwidth = 4.050000e-02; TDILevel = 24; END_GROUP = BAND_C
BEGIN_GROUP = BAND_B
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 1.160760e-02; effectiveBandwidth = 5.400000e-02; TDILevel = 21; END_GROUP = BAND_B
BEGIN_GROUP = BAND_G
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 8.856680e-03; effectiveBandwidth = 6.180000e-02; TDILevel = 21; END_GROUP = BAND_G
BEGIN_GROUP = BAND_Y
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 6.843070e-03; effectiveBandwidth = 3.810000e-02; TDILevel = 24; END_GROUP = BAND_Y
BEGIN_GROUP = BAND_R
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 1.020000e-02; effectiveBandwidth = 5.850000e-02; TDILevel = 14; END_GROUP = BAND_R
BEGIN_GROUP = BAND_RE
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 6.219960e-03; effectiveBandwidth = 3.870000e-02; TDILevel = 24; END_GROUP = BAND_RE BEGIN_GROUP = BAND_N
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 1.179710e-02; effectiveBandwidth = 1.004000e-01; TDILevel = 14; END_GROUP = BAND_N
BEGIN_GROUP = BAND_N2
ULLon = -120.89338929; ULLat = 47.79786951; ULHAE = 1241.94; URLon = -120.70434797; URLat = 47.79312915; URHAE = 648.78; LRLon = -120.70392538; LRLat = 47.73691772; LRHAE = 989.68; LLLon = -120.89205130; LLLat = 47.74075361; LLHAE = 1813.30; absCalFactor = 1.063780e-02; effectiveBandwidth = 8.890000e-02; TDILevel = 24; END_GROUP = BAND_N2 outputFormat = "GeoTIFF";
BEGIN_GROUP = IMAGE_1
satId = "WV03"; mode = "FullSwath"; scanDirection = "Forward"; CatId = "1040010021D24C00"; TLCTime = 2016-09-22T19:19:19.004450Z; numTLC = 2; TLCList = ( (0, 0.000000), (4336, 0.867200) firstLineTime = 2016-09-22T19:19:19.004450Z; avgLineRate = 5000.00; exposureDuration = 0.00020000; minCollectedRowGSD = 1.414; maxCollectedRowGSD = 1.415; meanCollectedRowGSD = 1.414; minCollectedColGSD = 1.345; maxCollectedColGSD = 1.346; meanCollectedColGSD = 1.346; meanCollectedGSD = 1.380; meanProductRowGSD = 1.415; meanProductColGSD = 1.341; meanProductGSD = 1.378; rowUncertainty = 29.99; colUncertainty = 11.97; minSunAz = 167.8; maxSunAz = 167.8; meanSunAz = 167.8; minSunEl = 41.7; maxSunEl = 41.7; meanSunEl = 41.7; minSatAz = 159.0; maxSatAz = 159.8; meanSatAz = 159.5; minSatEl = 67.9; maxSatEl = 68.4; meanSatEl = 68.1; minInTrackViewAngle = -16.7; maxInTrackViewAngle = -16.7; meanInTrackViewAngle = -16.7; minCrossTrackViewAngle = 10.5; maxCrossTrackViewAngle = 10.6; meanCrossTrackViewAngle = 10.5; minOffNadirViewAngle = 19.6; maxOffNadirViewAngle = 19.7; meanOffNadirViewAngle = 19.7; PNIIRS = 3.5; cloudCover = 0.003; resamplingKernel = "CC"; positionKnowledgeSrc = "R"; attitudeKnowledgeSrc = "R"; revNumber = 11444; END_GROUP = IMAGE_1 END;

  • Why can't you define the extent with extent(g)=extent(raster(xmn=...xmx=...,ymn=...,ymx=...)) with the correct extent boundaries? – Spacedman Dec 14 '17 at 22:30
  • What does command-line gdalinfo tell you about the tif? Are there any other metadata files with it? – Spacedman Dec 14 '17 at 22:31
  • I can't run gdalinfo becasue I do not have GDAL installed and cannot install it on this work computer. When I try to add the extents it gives this error: unable to find an inherited method for function ‘extent’ for signature ‘"missing"’ – steelhead Dec 14 '17 at 22:44
  • There are several associated files that arcmap might be using, including: .RPB, .IMD, .GEO, .EPH, and .ATT – steelhead Dec 14 '17 at 22:47
  • r=raster() ; extent(r)=extent(raster(xmn=1,xmx=2,ymn=1,ymx=2)) should work to reassign an extent. Note the construction of a raster within the extent call. Otherwise you get that "missing" error. – Spacedman Dec 14 '17 at 22:48
1

This suggests that this is not a good geotiff file. But you can fix it, as pointed out by Spacedman, with extent. I would do:

g <- raster('16SEP22191919-M1BS-057310833010_01_P001.tif' )
extent(g) <- extent(-120.893178348, -120.703696534, 47.7366916019, 47.797656863)

(Note that you cannot name the arguments), unless you would do

 extent(g) <- extent(x=c(-120.893178348, -120.703696534, 47.7366916019, 47.797656863))

Also, you are not getting an error, you are getting a warning. It should still plot.

I think you do have GDAL installed via rgdal. For GDALinfo you can do:

library(rgdal)
GDALinfo('16SEP22191919-M1BS-057310833010_01_P001.tif' )

Edit

Reading the imd file is something gdal should do, but apparently does not, at least not in current rgdal. To work around it, this works for the file you supplied:

x <- readLines('test.imd')
y <- data.frame(strsplit(unlist(strsplit(x[grep("ULLon", x)[1]], ";")), "="), stringsAsFactors=FALSE)
colnames(y) <- trim(y[1,])
y <- y[-1,]
ext <- extent(as.numeric(c(y$ULLon, y$LRLon, y$LRLat, y$ULLat)))

extent(g) <- ext
  • GDALinfo did work. This is what my Tiff shows:rows 4336 columns 10651 bands 8 lower left origin.x 0 lower left origin.y 0 res.x 1 res.y 1 ysign -1 oblique.x 0 oblique.y 0 driver GTiff projection +proj=longlat +datum=WGS84 +no_defs – steelhead Dec 14 '17 at 23:06
  • right, the file is not georeferenced – Robert Hijmans Dec 14 '17 at 23:08
  • So, is ArcMap using a different metafile to correctly import it? – steelhead Dec 14 '17 at 23:10
  • 1
    Perhaps there is a ".tfw" file ("world file") that ArcGIS uses? – Robert Hijmans Dec 14 '17 at 23:11
  • I suspect the .imd file might be the culprit (desktop.arcgis.com/en/arcmap/latest/manage-data/…) - can you view its contents and maybe edit them into your question? – Spacedman Dec 15 '17 at 8:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.