0

I am trying to create a raster with Sentinel 3 data. I found out that long/lat coordinates were in a different file than the data (chla). I am using this code to create one raster that I would like to use to extract chla at specific locations.

library(raster)
#> Loading required package: sp

longitude <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/geo_coordinates.nc",
    varname = "longitude"
  )
#> Loading required namespace: ncdf4
latitude <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/geo_coordinates.nc",
    varname = "latitude"
  )
chla <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/chl_oc4me.nc"
  )
#> Warning in .varName(nc, varname, warn = warn): varname used is: CHL_OC4ME
#> If that is not correct, you can set it to one of: CHL_OC4ME, CHL_OC4ME_err

extent(chla) <-
  c(min(values(longitude)), max(values(longitude)), min(values(latitude)), max(values(latitude)))

projection(chla) <-
  CRS("+proj=longlat +datum=WGS84")

wm <- rworldmap::getMap(resolution = "high")
wm <- crop(wm, extent(-180, 180, 45, 90))
#> Loading required namespace: rgeos
#> Warning in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func): spgeom1 and
#> spgeom2 have different proj4 strings
#> Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
#> unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
#> strings
wm <- spTransform(wm, CRSobj = CRS(proj4string(chla)))

I get weird results when plotting data as is:

plot(chla)
plot(wm, add = TRUE)

If I flip the data it looks better, but it is not quite there. The data is shifted toward the upper right.

plot(flip(chla, 2))
plot(wm, add = TRUE)

Anyone can spot the mistakes I make?

Update #2

This is what it should look like:

Example

This is the scatterplot between longitude and latitude.

enter image description here

Update #3

I think I managed to get it to work using the rasterize() function. However, I am still not sure if I do it the right way and how to adequately choose nrow and ncol in the rasterize() function.

library(raster)
#> Loading required package: sp

longitude <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/geo_coordinates.nc",
    varname = "longitude"
  )
#> Loading required namespace: ncdf4
latitude <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/geo_coordinates.nc",
    varname = "latitude"
  )
chla <-
  raster(
    "/media/data4tb/sentinel3/Data to Philippe/S3A_OL_2_WFR____20160608T142332_20160608T142532_20171101T060408_0119_005_096______MR1_R_NT_002.SEN3/chl_oc4me.nc"
  )
#> Warning in .varName(nc, varname, warn = warn): varname used is: CHL_OC4ME
#> If that is not correct, you can set it to one of: CHL_OC4ME, CHL_OC4ME_err

e <- extent(cbind(values(longitude), values(latitude)))
r <- raster(e, ncol = 1000, nrow = 1000) # How to choose ncol and nrow correctly?

# Remove NA values from chla, faster for the rasterize function.
xyz <- cbind(values(longitude), values(latitude), values(chla))
xyz <- na.omit(xyz)

x <- rasterize(xyz[, 1:2], r, xyz[, 3], fun = mean)
crs(x) <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs "

# World map

wm <- rworldmap::getMap(resolution = "high")
wm <- crop(wm, extent(-180, 180, 45, 90))
#> Loading required namespace: rgeos
#> Warning in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func): spgeom1 and
#> spgeom2 have different proj4 strings
#> Warning in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td,
#> unaryUnion_if_byid_false, : spgeom1 and spgeom2 have different proj4
#> strings
# wm <- spTransform(wm, CRSobj = CRS(proj4string(chla)))

plot(x)
plot(wm, add = TRUE)

5
  • Any chance of a link to the data? And what should it look like? Should all the data (green-orange bits of the raster) be over land or sea or something? I can't see why you think the second one looks better!
    – Spacedman
    Mar 1, 2019 at 7:54
  • I suspect that the reason the lat and long are supplied as a raster is because the data are not on a regular lat-long grid. To confirm this do a scatterplot of the values() of the longitude raster vs the values of the latitude raster - if that's not a regular grid then you cant force the data to a grid based on the min-max as you do now.
    – Spacedman
    Mar 1, 2019 at 8:16
  • Here a link for the data: file.io/zat2fL. I should not have color over the land. If I plot the data in Snap I have something different (I will update my post with what is expected). Mar 1, 2019 at 12:39
  • Correct url for the data: file.town/download/vhx1kdbec7t94dose4bml46t2 Mar 1, 2019 at 12:52
  • 1
    I think now you understand the data the question has shifted. How you manipulate the data will depend on your application. For example if you have standardised your analysis on a particular grid then you'd rasterise to that, if you have point data you might do nearest-neighbour to each of the points, or an average...
    – Spacedman
    Mar 1, 2019 at 17:00

1 Answer 1

1

In case people are still looking at this RE: "how to adequately choose nrow and ncol in the rasterize() function"

In order to set nrow and ncol for your "blank" raster 'r', you can look to the nrow and ncol of one of the 3 rasters you loaded in. This will then lead to your final raster product having the correct resolution etc. Note that - in my experience - these nrow and ncol values will differ per scene, so if you are batch processing in a loop for example, you will want to iteratively extract these values for each scene.

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.