I am working with R statistics. And I am trying to extract pixel values for each point in each Landsat bands, however, I am just getting NAs instead of getting pixel values

Does anyone have any other ideas about what I am doing wrong?

##point files
trainData_crs<- read.csv(file.csv)

##images loading and convert it to bricks
all_landsat_bands<-list.files ("file.tif")
img <- brick(all_landsat_bands)
##convert to same projection 

## conver~t csv to shp

coordinates(trainData_crs)<-~Latitude +Longitude
proj4string(trainData_crs) <- ("+proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
+no_defs ")

##extract pixel values 
trainData_bands<-raster::extract(x=img,y=trainData_sf, df=T)

I don't know how to supply my data, I don't see anywhere to attach it, but I am also new to posting here but I will try to give you the printout of the data.

##point files
trainData_crs<- read.csv(file.csv)                                                                       > class(trainData_crs)
[1] "data.frame"

> summary(trainData_crs)
       PlotType      Class          Latitude       Longitude    
 closedforest:22   Min.   :1.000   Min.   :17.81   Min.   :24.72  
 grass       :34   1st Qu.:2.000   1st Qu.:17.86   1st Qu.:24.98  
 lowforest   :24   Median :3.000   Median :17.92   Median :25.02  
 mediumforest:26   Mean   :2.904   Mean   :17.95   Mean   :25.02  
 shrubs      : 4   3rd Qu.:4.000   3rd Qu.:18.01   3rd Qu.:25.10  
 unknown     : 3   Max.   :7.000   Max.   :18.46   Max.   :25.16  
 water       : 1 

##images load and convert it to bricks
all_landsat_bands<-list.files ("file.tif")
img <- brick(all_landsat_bands)
> img <- brick(LandsatAprJune_2019_st)
> class(img)
[1] "RasterBrick"
[1] "raster"
> names(img)
[1] "B1"  "B2"  "B3"  "B4"  "B5"  "B6"  "B7"  "B10" "B11"

##convert to same projection 

coordinates(trainData_crs)<-~Latitude +Longitude
trainData_crs <- SpatialPoints(cbind(trainData_crs$Latitude, trainData_crs$Longitude), 
                            proj4string=CRS("+proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
+no_defs "))
trainData_crs_sp <- spTransform(trainData_crs, crs(img))

##confirm if points and image are in same coordinates
> crs(trainData_crs_sp)
CRS arguments:
 +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
> crs(img)
CRS arguments:
 +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m

##summary on transformed points
> summary(trainData_crs_sp)
Object of class SpatialPoints
                 min        max
coords.x1   177676.9   244145.4
coords.x2 12738109.3 12786443.7
Is projected: TRUE 
proj4string :
[+proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m
Number of points: 114

##extract pixel values 
trainData_bands<-raster::extract(img, trainData_crs_sp)

I have edited the code above to try using spTransform on the points but I am still getting NAs

Read data files

trainData_crs<- read.csv(".csv")

coordinates(trainData_crs)<-~Latitude +Longitude

set lat-long coordinate numbers

proj4string(trainData_crs) = CRS("+init=epsg:4326")

crs(trainData_crs) CRS arguments: +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

ransform them to the CRS of the image

trainData_crs_sp <- spTransform(trainData_crs, crs(img))

plot images

plotRGB(img, r = 4, g = 3, b = 2, stretch = "lin", scale=1000, axes = TRUE, main = "RGB composite image\n Landsat Bands 7, 3, 2") box(col = "black"); plot(trainData_crs_sp, add=T)


crs(img) crs(trainData_crs_sp)

crs(trainData_crs_sp) CRS arguments: +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m +no_defs crs(img) CRS arguments: +proj=utm +zone=34 +south +a=6378249.145 +b=6356514.966398753 +towgs84=-138,-105,-289,0,0,0,0 +units=m +no_defs

They seem to be in the same coordinate but when I plot them , they are not appearing in the same spot/they are not overlapping. I think this is why I am getting NAs because there is no overlap. any idea of the possibilities Why they are not overlapping if they are in the same coordinates now?

  • Hard to tell without seeing your data. But if the lat-long numbers in the CSV are really lat-long (in degrees) then setting the projection string to that UTM zone 34 CRS will mess things up bigly because they aren't UTM 34 numbers. But maybe they are. Can't tell without some info about your data. If you cant supply your data please at least edit and show some summary info about the data objects you are creating, eg plots or prints. – Spacedman Nov 7 '19 at 18:53
  • It looks like you may be just assigning a coordinate system and not actually transforming the data. After coercing to an sp object, try using spTransform on your points. – Jeffrey Evans Nov 7 '19 at 18:56

The numbers in your point CSV look like degrees in lat-long. This is (999 times out of a 1000) EPSG coordinate system 4326. After reading the data in you need to assign that CRS to it before you transform it to anything else. Note there is a great difference between assigning a coordinate system to some data (which does not change the coordinate numbers) and transforming it to another system (which does change the numbers).

trainData_crs<- read.csv(file.csv)

## convert csv to shp

coordinates(trainData_crs)<-~Latitude +Longitude

at this point you've got a spatial data set with those lat-long coordinate numbers. Tell R these are in lat-long:

proj4string(trainData_crs) = CRS("+init=epsg:4326")

Now you should be able to transform them to the CRS of the image:

trainData_crs_sp <- spTransform(trainData_crs, crs(img))

Plot the image [plot(img)] and add the points [plot(trainData_crs_sp, add=TRUE)] to check the points overlap.

| improve this answer | |

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.