[R-sig-Geo] Problem in extracting data from GRIB files

[Michael Sumner]

On Sat, 29 Jul 2017 at 02:21 Miluji Sb <milujisb at gmail.com> wrote:

> Dear all,
>
> I have a set of coordinates for cities:
>
> structure(list(lon = c(3.7174243, 3.2246995, 2.928656, 33.3822764,
> 10.40237, 24.7535746, 7.2619532, -0.370679, -4.486076, -4.097899
> ), lat = c(51.0543422, 51.209348, 51.21543, 35.1855659, 55.403756,
> 59.4369608, 43.7101728, 49.182863, 48.390394, 47.997542)), .Names =
> c("lon",
> "lat"), na.action = structure(c(5L, 6L, 31L, 55L, 59L, 61L, 67L,
> 68L), .Names = c("5", "6", "31", "55", "59", "61", "67", "68"
> ), class = "omit"), row.names = c(1L, 2L, 3L, 4L, 7L, 8L, 9L,
> 10L, 11L, 12L), class = "data.frame")
>
> I am trying to extract climatic data from GLDAS climatic data (1 degree x 1
> degree) GRIB files with the following code:
>
> # Convert to spatial points
> pts_dams <- SpatialPoints(lonlat,proj4string=CRS("+proj=longlat"))
>
> # Extract
>  filelist<-list.files(pattern=".grb$")
>
>
This part is a bit wrong:


>     data<-stack(filelist[1])
>     data at layers<-sapply(filelist, function(x) raster(x,band=12))
>

Better would be

rdata <-  raster::stack(lapply(filelist, function(x) raster(x,band=12)))

For the rest we will need to see details about the data, so please do

print(rdata)

and send the resulting print-out.

 A minor thing "data" is already a function, so it's advisable not to use
it as a name - if  only to reduce possible confusion.

Also, it's not obvious that the process by which raster(file) goes though
(it passes down to rgdal::readGDAL) will result in a correct interpretation
as the data source was intended, since GRIB is a domain-specific format for
time-series and volume-slice series data and the required translation to
the GDAL and raster data models is not always straight-forward.

It totally can work though, I do it routine with many sources but they all
need some oversight to ensure that upfront translation is complete and
appropriate.

It may just require setting raster::extent and/or raster::projection, but
there's no way to know without exploration.

Cheers, Mike.



>     data_dams_size<-extract(data,pts_dams)
>     DF_data_dams_size<- as.data.frame(data_dams_size)
>
> Unfortunately, I get mostly NAs for the data, it seems that there's an
> issue with the CRS projections for the city coordinates. Is there a
> specific projection for city level coordinates? Or am I doing something
> completely wrong? Thank you!
>
> Sincerely,
>
> Milu
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
-- 
Dr. Michael Sumner
Software and Database Engineer
Australian Antarctic Division
203 Channel Highway
Kingston Tasmania 7050 Australia

	[[alternative HTML version deleted]]