CollinWoo / daynight-Q10

MIT License
2 stars 0 forks source link

Reading grib data #24

Open CollinWoo opened 3 years ago

CollinWoo commented 3 years ago

I don't know what the data in the .grib files are trying to convey. The projection type seems to be stereographic (https://proj.org/operations/projections/stere.html), so I have no idea if it extracts data based on latitude and longitude. My function for extracting moisture data from specific coordinates seems to only work for very specific values and returns NA for all others. Having a better understanding of stereographic projections and how to convert them to other formats (namely, latlon) would be a good start.

Here is an image of the information in the raster object extracted from one of the .grib files: image

bpbond commented 3 years ago
library(rgdal)
x <- readGDAL("~/Data/NASA-USDA_SoilMoisture/L03/20100110_20100112.as1.grib")
plot(x)

Screen Shot 2020-10-09 at 6 06 59 AM

That seems good.

bpbond commented 3 years ago

I played around a bit more:

library(raster)
y <- brick("~/Data/NASA-USDA_SoilMoisture/L03/20100110_20100112.as1.grib")
xy <- cbind(50, seq(-80, 80, by=20))  # crosses the Sahara so should be easy to see if correct
sp <- SpatialPoints(xy, proj4string = CRS("+proj=longlat +datum=WGS84"))
extract(y, sp)
     X20100110_20100112.as1
 [1,]                     NA
 [2,]                     NA
 [3,]                     NA
 [4,]                     NA
 [5,]                     NA
 [6,]                     NA
 [7,]                     NA
 [8,]                     NA
 [9,]                     NA
Warning message:
In .local(x, y, ...) : Transforming SpatialPoints to the CRS of the Raster

Hmm well THAT didn't work. I don't have time to pursue this further right now, but leaving this note for when I (or one of you) try again.

bpbond commented 3 years ago

OK this is really frustrating. @jinshijian do you think Dalei would have advice on how to do this?

CollinWoo commented 3 years ago

I think it may have something to do with the projection since there is a warning message that prints beforehand. This stack overflow post may be relevant but I'm not sure how to implement it: https://stackoverflow.com/questions/63374488/since-update-of-sp-package-i-get-a-warning-by-calling-a-spcrs-definition.