2

I want to plot sea surface temperature data from a regular grid but can't find the proper way to do it. My data comes in nc format and can be downloaded from http://www.nodc.noaa.gov/SatelliteData/pathfinder4km/

I use this R code to access data but the problem appears when trying to plot

library("ncdf")
download.file("ftp://ftp.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.2/2003/20030715000715-NODC-L3C_GHRSST-SSTskin-AVHRR_Pathfinder-PFV5.2_NOAA17_G_2003196_night-v02.0-fv02.0.nc", destfile="sst.nc")
data=open.ncdf("sst.nc")
x <- get.var.ncdf(data,"lon")
y <- get.var.ncdf(data,"lat")
sst=get.var.ncdf(data,"sea_surface_temperature")

filled.contour(x,y,sst, color = terrain.colors, asp = 1)

And then get this error message

Error en filled.contour(x, y, sst, color = terrain.colors, asp = 1) : increasing 'x' and 'y' values expected

I think the problem comes from y coordinate, latitude runs from 90 to -90. I've seen some questions on stackoverflow in creating a new grid with akima package but it should not be necessary in this case.

Here you can find a summary of data file

http://ubuntuone.com/1mIdYVqoePn24gKQbtXy7K

Thanks in advance for your help

SOLVED

Thanks to Paul Hiemstra

The point was not read lat-lon values from data set but to know the i,j coordinates of data points in the matrix and then select the geographic area I want to plot. Following are the commands that work for me:

library("ncdf")
download.file("ftp://ftp.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.2/2003/20030715000715-NODC-L3C_GHRSST-SSTskin-AVHRR_Pathfinder-PFV5.2_NOAA17_G_2003196_night-v02.0-fv02.0.nc", destfile="sst.nc")
data=open.ncdf("sst.nc")
sst=get.var.ncdf(data,"sea_surface_temperature")
x = seq(1, 8640, length.out = nrow(sst))         # Matrix dimension 8640x4320
y = seq(1, 4320, length.out = ncol(sst))

sst1 <- sst[c(1000:1500),c(1000:1500)]           # Subsetting a region
x = seq(1, 500, length.out = nrow(sst1))
y = seq(1, 500, length.out = ncol(sst1))

png(filename="sst.png",width=800,height=600,bg="white")
filled.contour(x,y,sst1, color = terrain.colors, asp = 1)
dev.off()

Now I have to figure out how to label the plot with longtiude and latitude at x-y coordinates.

pacomet
  • 5,011
  • 12
  • 59
  • 111

2 Answers2

3

Although my question was solved with the help of Paul Hiemstra I've still investigated on plotting my netcdf data and found another thread in Stackoverflow that has also helped me. It uses image instead of filled.contour.

You can find the thread at The variable from a netcdf file comes out flipped

Now, this is the code I am using to plot SST data:

library("ncdf")
download.file("ftp://ftp.nodc.noaa.gov/pub/data.nodc/pathfinder/Version5.2/2000/20000107010122-NODC-L3C_GHRSST-SSTskin-AVHRR_Pathfinder-PFV5.2_NOAA14_G_2000007_night-v02.0-fv01.0.nc", destfile="sst.nc")
data=open.ncdf("sst.nc")


sst=get.var.ncdf(data,"sea_surface_temperature")
lon=get.var.ncdf(data,"lon")
lat=get.var.ncdf(data,"lat")

data$dim$lon$vals -> lon
data$dim$lat$vals -> lat
lat <- rev(lat)
sst <- sst[,ncol(sst):1]
png(filename="sst2.png",width=1215,height=607,bg="white")
image(lon, lat, sst, zlim=c(270,320), col = heat.colors(37))
library("sp", lib.loc="/usr/lib/R/site-library")
library("maptools", lib.loc="/usr/lib/R/site-library")
data(wrld_simpl)
plot(wrld_simpl, add = TRUE)
dev.off()

that leads to this image: enter image description here

Community
  • 1
  • 1
pacomet
  • 5,011
  • 12
  • 59
  • 111
2

The problem is probably the following. Your x and y variables are of the same size of sst, i.e. they provide maps of x and y coordinates. What the filled_contour function needs is not these maps, but more simply the x and y coordinates associated with the rows and columns in sst with increasing x values.

Paul Hiemstra
  • 59,984
  • 12
  • 142
  • 149
  • Thanks for answering Paul. So, what should I provide to filled_contour or image.plot? What if I want to show latitude and longitude in axis labels? – pacomet May 03 '13 at 09:21
  • You should provide `filled.contour` with two vectors mapping the rows and columns to latitude and longitude. You can even use `image` without specifying latitude and longitude. – Paul Hiemstra May 03 '13 at 09:23
  • Hi again, I'm not an experienced R user and so I'm not sure how can I map rows and columns as you say. Can you point me to any info I can read to look for that? Thanks – pacomet May 03 '13 at 12:00
  • The documentation of `filled.contour` specifies exactly what you need: `?filled.contour`. – Paul Hiemstra May 03 '13 at 12:01
  • Thanks for your help Paul, it's now working. Now I have to figure out how to use lat-lon to label axis. – pacomet May 06 '13 at 09:43