How to find the 4 coordinates (lat-long) of a pixel (cell) in a global file?

Question

I am using a climate variable that can be downloaded from here:

  ftp://sidads.colorado.edu/pub/DATASETS/nsidc0301_amsre_ease_grid_tbs/global/ 

This file is a binary (matrix) file with 586 lines and 1383 columns (global map). I would like to know the 4 coordinates (lat-long) of a pixel (cell)`.

more info about the file :

These data are provided in EASE-Grid projections global cylindricalat 25 km resolution, are two-
      byte 
    Spatial Coordinates:
    N: 90°     S: -90°     E: 180°     W: -180°    

use the raster package and convert the data to raster objects:

 file<- readBin("ID2r1-AMSRE-ML2010001A.v03.06H", integer(), size=2,  n=586*1383, signed=T)
 m = matrix(data=file,ncol=1383,nrow=586,byrow=TRUE)
 r = raster(m, xmn=-180, xmx=180, ymn=-90, ymx=90)

Now the file is a properly spatially referenced object, but without a full specification of the cylindrical projection used you can't get back to lat-long coordinates.

There some more info here http://nsidc.org/data/ease/tools.html including a link to some grids that have the lat-long of grid cells for that grid system:

ftp://sidads.colorado.edu/pub/tools/easegrid/lowres_latlon/

      MLLATLSB  "latitude"
      MLLonLSB  "longitude"

so for example we can create a raster of latitude for the cells in your data grid:

> lat <- readBin("MLLATLSB",integer(), size=4,  n=586*1383, endian="little")/100000
> latm = matrix(data=lat,ncol=1383,nrow=586,byrow=TRUE)
> latr = raster(latm, xmn=-180, xmx=180, ymn=-90, ymx=90)

and then latr[450,123] is the latitude of cell [450,123] in my data. Repeat with MLLONLSB for longitude.

However this is not sufficient(one lat and long for each pixel) as I would like to compare with ground based measurements so I need to define my region which correspond to this pixel (25 * 25 km or 0.25 degrees). for this purpose I have to know the 4 coordinates (lat-long) of that pixel (cell). Thanks for any help

Answer 1

EASE GRID uses a Global Cylindrical Equal-Area Projection defined by EPSG 3410, which is metric. As long as I see it, spatial extent should be provided in meters, not geographic coordinates. From here we see that map extent coordinates are:

xmin: -17609785.303313

ymin: -7389030.516717

xmax: 17698276.686747

ymax: 7300539.133283

So slightly changing your code we have this

library(raster)
library('rgdal')
wdata <- 'D:/Programacao/R/Raster/Stackoverflow'
wshp <- 'S:/Vetor/Administrativo/Portugal'
#setwd(wdata)
file <- readBin(file.path(wdata, "ID2r1-AMSRE-ML2010001D.v03.06H"),
               integer(), size=2,  n=586 * 1383, signed=T)
m <- matrix(data = file, ncol = 1383, nrow = 586, byrow = TRUE)
-17609785.303313 -7389030.516717 17698276.686747 7300539.133283
rm <- raster(m, xmn = -17609785.303313, xmx = 17698276.686747,
             ymn = -7389030.516717, ymx = 7300539.133283)
proj4string(rm) <- CRS('+init=epsg:3410')

> rm
class       : RasterLayer 
dimensions  : 586, 1383, 810438  (nrow, ncol, ncell)
resolution  : 25530.05, 25067.53  (x, y)
extent      : -17609785, 17698277, -7389031, 7300539  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:3410 +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +a=6371228 +b=6371228 +units=m +no_defs 
data source : in memory
names       : layer 
values      : 0, 3194  (min, max)

writeRaster(rm, file.path('S:/Temporarios', 'easegrridtest.tif'), overwrite = TRUE)
plot(rm, asp = 1)

We can now overlay some spatial data

afr <- readOGR(dsn = file.path(wdata), layer = 'Africa_final1_dd84')
proj4string(afr) <- CRS('+init=epsg:4326') # Asign projection 
afr1 <- spTransform(afr, CRS(proj4string(rm)))

plot(afr1, add = T)

Now you can start playing with extract extent of your ROI, possibly with extent()

I'm not happy with the spatial adjustment. With this approach I got a huge error. I'm not sure about the position error but it is described for this product. Maybe something with the extent parameters.

Since you are interested in use it against ground measures, you could polygonize your ROI and use it in your GPS or GIS.

Also you could get extent of cells of interest with something like:

Choose an aproximate coordinate, identify cell and get the extent:

cell <- cellFromXY(rm, matrix(c('x'= -150000, 'y' =200000), nrow = 1, byrow = T))
r2 <- rasterFromCells(rm, cell, values=TRUE)
extent(r2)
class       : Extent 
xmin        : -172759.8 
xmax        : -147229.7 
ymin        : 181362 
ymax        : 206429.6 

And maybe identify a ROI (single cell) in a map (or plot)

cell <- cellFromXY(rm, matrix(c('x'= -1538000, 'y' =1748000), nrow = 1, byrow = T))
r2 <- rasterFromCells(rm, cell, values=TRUE)
r2p <- as(r2, 'SpatialPolygons')
extr2 <- extent(r2) + 300000
plot(rm,  col = heat.colors(6), axes = T, ext = extr2)
plot(afr1, add = T, col = 'grey70')
plot(r2p, add = T)

For a particular area

And assuming that by cell you mean column and row values, you could proceed with raster::cellFromRowCol

cell2 <- cellFromRowCol(rm, rownr = mycellnrow, colnr = mycellnrow)
r3 <- rasterFromCells(rm, cell2, values=TRUE)
r3p <- as(r3, 'SpatialPolygons')
extr3 <- extent(r3) + 3000000

In this particular case, 123 and 450 seems to be far from any continental area...

Hope it helps.

More info on AMSR-E/Aqua Daily Gridded Brightness Temperatures here