1

I have an obstacle which I thus far wasn't able to overcome. I have gigabytes of Sentinel 2 imagery from which I need to extract pixel images. Firstly, I have converted the images to .tif files. After that, in some cases, the extract function works as I expect it to:

        image=raster("E:/S2A_USER_MSI_L2A_TL_SGS__20160317T153554_ 
        A003837_T33TWJ_B8A_20m.tif")
        typeof(image)

    [1] "S4

        str(image)

            Formal class 'RasterLayer' [package "raster"] with 12 slots
              ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
              .. .. ..@ name        : chr "E:\\0103_2404\\S2A_USER_PRD_MSIL2A_PDMC_20160309T094150_R136_V20160308T093139_20160308T093139.SAFE\\GRANULE\\S2"|
        __truncated__
              .. .. ..@ datanotation: chr "INT2U"
              .. .. ..@ byteorder   : chr "little"
              .. .. ..@ nodatavalue : num -Inf
              .. .. ..@ NAchanged   : logi FALSE
              .. .. ..@ nbands      : int 1
              .. .. ..@ bandorder   : chr "BIL"
              .. .. ..@ offset      : int 0
              .. .. ..@ toptobottom : logi TRUE
              .. .. ..@ blockrows   : int 1
              .. .. ..@ blockcols   : int 5490
              .. .. ..@ driver      : chr "gdal"
              .. .. ..@ open        : logi FALSE
              ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
              .. .. ..@ values    : logi(0) 
              .. .. ..@ offset    : num 0
              .. .. ..@ gain      : num 1
              .. .. ..@ inmemory  : logi FALSE
              .. .. ..@ fromdisk  : logi TRUE
              .. .. ..@ isfactor  : logi FALSE
              .. .. ..@ attributes: list()
              .. .. ..@ haveminmax: logi TRUE
              .. .. ..@ min       : num 0
              .. .. ..@ max       : num 65535
              .. .. ..@ band      : int 1
              .. .. ..@ unit      : chr ""
              .. .. ..@ names     : chr "S2A_USER_MSI_L2A_TL_SGS__20160317T153554_A003837_T33TWJ_B8A_20m"
              ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
              .. .. ..@ type      : chr(0) 
              .. .. ..@ values    : logi(0) 
              .. .. ..@ color     : logi(0) 
              .. .. ..@ names     : logi(0) 
              .. .. ..@ colortable: logi(0) 
              ..@ title   : chr(0) 
              ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
              .. .. ..@ xmin: num 5e+05
              .. .. ..@ xmax: num 609780
              .. .. ..@ ymin: num 4790220
              .. .. ..@ ymax: num 4900020
              ..@ rotated : logi FALSE
              ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
              .. .. ..@ geotrans: num(0) 
              .. .. ..@ transfun:function ()  
              ..@ ncols   : int 5490
              ..@ nrows   : int 5490
              ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
              .. .. ..@ projargs: chr "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
              ..@ history : list()
              ..@ z       : list()

        pixel_values=(na.omit(extract(image,azo_plohe, cellnumbers=TRUE, 
        weights=TRUE, df=TRUE)))


     ID   CELL pxl_value  Weight
1   320 5625881 0   0.034532374
2   320 5625882 0   0.082014388
3   320 5625883 0   0.027338129
4   320 5631371 0   0.135251799
5   320 5631372 0   0.143884892
6   320 5631373 0   0.122302158
7   320 5636861 0   0.125179856
8   320 5636862 0   0.143884892
9   320 5636863 0   0.110791367

In this case my first 9 pixels are on the areas where it is all black so the values are zero, but the function works. However, in some other cases when I use this same function on another image:

image=raster("E:/SS2A_USER_MSI_L2A_TL_SGS__ 20160308T165710_A003708_T34TCN_B12_20.
typeof(image)

  [1] "S4

            str(image)

Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "E:\\0103_2404\\S2A_USER_PRD_MSIL2A_PDMC_20160309T094150_R136_V20160308T093139_20160308T093139.SAFE\\GRANULE\\S2"| __truncated__
  .. .. ..@ datanotation: chr "INT2U"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ bandorder   : chr "BIL"
  .. .. ..@ offset      : int 0
  .. .. ..@ toptobottom : logi TRUE
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 5490
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
  .. .. ..@ values    : logi(0) 
  .. .. ..@ offset    : num 0
  .. .. ..@ gain      : num 1
  .. .. ..@ inmemory  : logi FALSE
  .. .. ..@ fromdisk  : logi TRUE
  .. .. ..@ isfactor  : logi FALSE
  .. .. ..@ attributes: list()
  .. .. ..@ haveminmax: logi TRUE
  .. .. ..@ min       : num 0
  .. .. ..@ max       : num 65535
  .. .. ..@ band      : int 1
  .. .. ..@ unit      : chr ""
  .. .. ..@ names     : chr "S2A_USER_MSI_L2A_TL_SGS__20160308T165710_A003708_T34TCN_B12_20m"
  ..@ legend  :Formal class '.RasterLegend' [package "raster"] with 5 slots
  .. .. ..@ type      : chr(0) 
  .. .. ..@ values    : logi(0) 
  .. .. ..@ color     : logi(0) 
  .. .. ..@ names     : logi(0) 
  .. .. ..@ colortable: logi(0) 
  ..@ title   : chr(0) 
  ..@ extent  :Formal class 'Extent' [package "raster"] with 4 slots
  .. .. ..@ xmin: num 3e+05
  .. .. ..@ xmax: num 409800
  .. .. ..@ ymin: num 4690200
  .. .. ..@ ymax: num 4800000
  ..@ rotated : logi FALSE
  ..@ rotation:Formal class '.Rotation' [package "raster"] with 2 slots
  .. .. ..@ geotrans: num(0) 
  .. .. ..@ transfun:function ()  
  ..@ ncols   : int 5490
  ..@ nrows   : int 5490
  ..@ crs     :Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=utm +zone=34 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
  ..@ history : list()
  ..@ z       : list()

           pixel_values=(na.omit(extract(image,azo_plohe, cellnumbers=TRUE, 
            weights=TRUE, df=TRUE)))

Error in names(x) <- value : 'names' attribute [4] must be the same length as the vector [0]

I do not understand why i can not execute extract function on this 2 seemengly same .tif files. So far i wasn't able to see any pattern why some images pixels are normally extracted and why some aren't. I know that if i remove the

df=TRUE

I do not get the error. But the thing is that i need this part for the rest of my code which can be found on this link (How to get feature values from multiple polygons of shapefile while using extract function on raster). IF I remove the "df=TRUE" I get another error

Error in .doXYFromCell(object@ncols, object@nrows, e@xmin, e@xmax, e@ymin,  : 
Not compatible with requested type: [type=NULL; target=double].

Sidenote, some images return null values because they are not in my area of interest so the extract works OK on SOME of them. But, I have noticed that images that give error when trying to extract are all outside of my area of interes. This is my newest clue regarding the problem of extraction.

  • Seems to be a problem with a long path name... But I'm not sure... Why convert S-2A images to tiff ? rgdal supports jpg-2000 (if you update gdal) – aldo_tapia Nov 7 '17 at 18:41
  • Hmmm, the paths are always long and should always be the same beacuse of the structure of folders when downloading Sentinel imagery (earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/… ) . If this were the problem then i shouldnt be able to extract pixels from any of the images right? I converted to tif because I couldn't coarse R into reading jp2 rasters with "raster(x)" function, even though i have several gdal jp2 drivers. – davor korman Nov 7 '17 at 18:47
  • Yes, I know S-2*.SAFE folder structure... But for many softwares is an issue (also, working in windows). Only in SNAP everything works fine. I'm an OS X user and this answer helped me to set-up gdal correctly in R to open S-2 images (S-2A and S-2B) – aldo_tapia Nov 7 '17 at 18:56
  • @davorkorman, I was also converting to GeoTIFF to work with S2 data in R. After uploading rgdal, I am able to use the jp2 format in R. Which version are you using? It may save you a step in your methodology. Anyways, did you try to save your tiff S2 files in a new path, out of the .SAFE structure? I usually do that and when extracting pixel values I do not face any problem. – GCGM Nov 8 '17 at 9:22
  • I am using 2.2.1. on Windows 10. Then it is possible that in future I might save a step with the conversion. I have found out that if i loose the "df=TRUE" in my code, I do not get the error "Error in names(x) <- value : 'names' attribute [4] must be the same length as the vector [0]"., but i get this error "Error in .doXYFromCell(object@ncols, object@nrows, e@xmin, e@xmax, e@ymin, : Not compatible with requested type: [type=NULL; target=double].". This error is a part of a bigger code I posted as another question. If you are interested in further helping I will post it bellow this coment. – davor korman Nov 8 '17 at 10:48

Your Answer

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

Browse other questions tagged or ask your own question.