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;

Improve this question
5
  • Why can't you define the extent with extent(g)=extent(raster(xmn=...xmx=...,ymn=...,ymx=...)) with the correct extent boundaries?
    – Spacedman
    Dec 14, 2017 at 22:30
  • What does command-line gdalinfo tell you about the tif? Are there any other metadata files with it?
    – Spacedman
    Dec 14, 2017 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, 2017 at 22:44
  • There are several associated files that arcmap might be using, including: .RPB, .IMD, .GEO, .EPH, and .ATT
    – steelhead
    Dec 14, 2017 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, 2017 at 22:48

1 Answer 1

Reset to default
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
Improve this answer
7
  • 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, 2017 at 23:06
  • right, the file is not georeferenced
    – Robert Hijmans
    Dec 14, 2017 at 23:08
  • So, is ArcMap using a different metafile to correctly import it?
    – steelhead
    Dec 14, 2017 at 23:10
  • 1
    Perhaps there is a ".tfw" file ("world file") that ArcGIS uses?
    – Robert Hijmans
    Dec 14, 2017 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, 2017 at 8:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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