Add long/lat points to map of landsat data

Question

I have a RasterBrick for a landsat image:

landsat_stack_brick
#class      : RasterBrick 
#dimensions : 7711, 7551, 58225761, 7  (nrow, ncol, ncell, nlayers)
#resolution : 30, 30  (x, y)
#extent     : 576885, 803415, 2281485, 2512815  (xmin, xmax, ymin, ymax)
#crs        : +proj=utm +zone=4 +datum=WGS84 +units=m +no_defs 
#source     : r_tmp_2022-11-19_144142_2384_62790.grd 
#names      :    B1,    B2,    B3,    B4,    B5,    B6,    B7 
#min values :     1,     1,     1,     4,  1780,  6138,  7110 
#max values : 60301, 61806, 61527, 62936, 62293, 65454, 65454

I want to add points in longitude/latitude to the plot of landsat_stack_brick

 library(raster)

 Long_lat_df  = na.omit(water_long_lat)
 coordinates(Long_lat_df ) <- ~ water_long + water_lat
 proj4string(Long_lat_df ) <- CRS("+init=epsg:4326")
 y <- spTransform(geo_points, crs(landsat_stack_brick))

 plot(landsat_stack_brick)
 points(y)

The points should be all along the coast lines of the islands, but this is how it looks when I try to add them to the map:

The code used thus far is: [EDIT: Applying corrections mentioned by John Polo]

water_chem <- read.csv("data.csv")

###Extract all first unique values obtained for latitude)
water_lat = unique(water_chem$Latitude)
###Extract all corresponding longitutude values for unique latitud ###values found
water_long= water_chem$Longitude[unique(water_chem$Latitude)]

x <- water_long_lat
x = na.omit(x)
coordinates(x) <- ~ water_long + water_lat
sbux_sf <- st_as_sf(x, coords = c("water_long", "water_lat"),  crs = landsat_stack@crs)

sbux_sd

structure(list(geometry = structure(list(structure(c(-156.63776, 
21.0133265), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, 
21.0134451), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, 
21.0135265), class = c("XY", "POINT", "sfg")), structure(c(-156.63776, [...]

[...], precision = 0, bbox = structure(c(xmin = -156.6407, ymin = 20.77606, 
xmax = -156.63776, ymax = 21.014317), class = "bbox"), crs = structure(list(
    input = "WGS 84", wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",6326]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433],\n        ID[\"EPSG\",8901]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]],\n    USAGE[\n        SCOPE[\"unknown\"],\n        AREA[\"World.\"],\n        BBOX[-90,-180,90,180]]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
124L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(integer(0), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = character(0)))

It would be really helpful to know what your `sbx_sf` looks like. Add the output of `dput(sbx_sf)` to your question. If that produces a lot of lines, you can use `dput(head(sbx_sf))` instead. Also, whenever you import libraries, you should include that with the code you share. It looks like you're at least using `sf` and `raster`, but maybe something else too? Add that also. And don't put pictures of data in your questions. Where you where you typed `landsat_stack_brick`, just paste that and the result in the code, not a picture of it. Pictures of maps that illustrate the problem are fine. — John Polo, Nov 20 '22 at 13:42
@John Polo, thanks for the suggestions. I added all the points you mentioned. — Pablo Rodriguez, Nov 20 '22 at 14:02
Thanks for making some edits, but I'm still confused with your workflow. Ideally, you should provide the code needed for a [minimum reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) and it should be in order and continuous. This code jumps around and I see sbux_sf and sbux_sd with no step in between, for example. And is the stuff after sbux_sd supposed to be the `dput` result? Something is wrong there. I say this because of the "[...]". That isn't part of a `dput` result. Use "#" in code for explanations if needed — John Polo, Nov 20 '22 at 14:24
I think one of your problems is the CRS you're using. You have EPSG:32604 in the Landsat image and EPSG:4326 in other places. I know you're trying to reproject, but if you don't have the right EPSG on the object to begin with, the reproject is not going to work. — John Polo, Nov 20 '22 at 14:30
Further thought tells me that unique of the lat is just going to give you one lat value from the source of your data and that's going to screw up a bunch of the vectors. Why are you doing that? If fixing that doesn't work you should probably provide the `dput` of the object that you create with `read.csv( water.csv)`, not sbx_df, since the water is the starting data for the polygon. — John Polo, Nov 20 '22 at 14:39

score 1 · Accepted Answer · answered Nov 21 '22 at 02:57

It is not clear what you have done exactly (please show one minimal, but complete, script, not bits and pieces). How did you create landsat_stack_brick?

This piece of your code is clearly wrong

water_lat = unique(water_chem$Latitude)
water_long= water_chem$Longitude[unique(water_chem$Latitude)]

You probably want

water <- unique(na.omit(water_chem[, c("Longitude", "Latitude")])

A complete workflow could like this (using "terra" as that has replaced the "raster" package):

library(terra)
r <- rast("landsat.tif")
water <- unique(na.omit(water_chem[, c(Longitude", "Latitude")])
w <- vect(water, c("Longitude", "Latitude"), crs="+proj=longlat")
w <- project(w, crs(r))

plotRGB(r, 3,4,5)
points(w)

Add long/lat points to map of landsat data

1 Answers1