wrf_raster wrong projection

[bakamotokatas]

Hello @Schuch666, wrf_raster works for me now. But,

wrf_raster is getting the projection of the wrf wrong (seems to reverse longitude and latitude). I think it's probably because my gdal library version is higher than 3.0. You can see the sample image below. This image actually needs to sit on Turkey.

I had created a similar issue for a similar GIS project. I think it had the same problem. They solved according to this link. But I couldn't make any changes in the code for a solution because I don't have much idea about it.

My example code;

library(raster)
library(ncdf4)
library(EmissV)
library(eixport)
library(leaflet)
library(htmlwidgets)
raster<-wrf_raster(file = "/home/bakamotokatas/Build_WRF/WRF-4.3-ARW/run/wrfout_d01_2020-02-16_00:00:00",name = "HGT")
plot(raster$HGT_2020.02.16_00.00.00)
raster<-projectRaster(raster, crs='+proj=longlat +datum=WGS84')
pal <- colorNumeric(c("forestgreen", "yellow", "red"), values(raster$HGT_2020.02.16_00.00.00),na.color = "transparent")
leaflet() %>%
  addTiles() %>%
  addRasterImage(raster$HGT_2020.02.16_00.00.00, colors = pal, opacity = 0.7) %>%
  addLegend(pal = pal, 
            values = values(raster$HGT_2020.02.16_00.00.00),
            title = "Height (m)")

[Schuch666]

Hi @bakamotokatas

Sorry to hear this new issue, I can't reproduce this problem but it doesn't seem to be related to projection but in the order of some arguments passed to the project function that is internally imported from the rgdal package. In next versions of the eixport this process will be performed using functions from sp/sf packages.

As a workaround, I added an argument that switches the lat for lon, you will need to set reverse = TRUE:

...
raster<-wrf_raster(file = "/home/bakamotokatas/Build_WRF/WRF-4.3-ARW/run/wrfout_d01_2020-02-16_00:00:00",name = "HGT",reverse = T) 
...

Let me know if it works for you

[bakamotokatas]

Thanks for the workaround @Schuch666. Actually, I tried swapping y with x too. But it still doesn't fit the map like in the picture.

I tried to play with this part, without reversing y and x.

from

  xmn <- projcoords[1,1] - dx/2  # Left border
  ymx <- projcoords[1,2] + dy/2   # upper border

to

  xmn <- projcoords[1, 2] - dx/2
  ymx <- projcoords[1, 1] + dy/2

it approached the area where it should be, but still further east and south of the place where it should be; as it is seen in the image.

I tried adding dx and dy with different combinations, but I still couldn't fit.

additionally, the code gives this warning while drawing.

Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition
3: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs
4: In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum World Geodetic System 1984 in Proj4 definition

[Schuch666]

Hi @bakamotokatas

After some testing I changed the wrf_raster function to use both rgdal and sf to do this transformation (rgdal will be removed in future versions), In the following example there is the use of this function with rgdal and using sf (the commented parts are only needed to plot the map).

library(eixport)
library(raster)
# library(maps)
# library(hackWRF)
file   = paste(system.file("extdata", package = "eixport"),"/wrfinput_d01", sep="")
# mapa   <- raster::shapefile(paste(system.file("extdata", package = "hackWRF"),"/BR_estates.shp",sep=""),verbose=F)
Rrgdal <- wrf_raster(file,'XLAT',latlon = T)
plot(Rrgdal, legend.shrink = 0.98,legend.width = 3, main = 'XLAT rgdal::project')
# map(mapa, add=TRUE, lwd=1, interior = FALSE, col = "black")

Rsf    <- wrf_raster(file,'XLAT',latlon = T,use_sf = T)
plot(Rsf, legend.shrink = 0.98,legend.width = 3, main = 'XLAT sf::st_transform')
# map(mapa, add=TRUE, lwd=1, interior = FALSE, col = "black")

these warnings related to projection are not a problem in general.

Let me know if this code is working and if this changes solves the problem, if it doesn't, I'll need more data to reproduce the problem.

[bakamotokatas]

Hi @Schuch666, really sorry for wasting your time. It's been my fault from the start. I didn't realize why the comma was there because I didn't examine it in detail. Now, when I examined it in detail, I realized it and you only need to correct the place below.

from

  if (length(dim(inNCLon)) == 3) {
    x <- as.vector(inNCLon[ncol(inNCLon):1])
    y <- as.vector(inNCLat[ncol(inNCLat):1])
  }

to

  if (length(dim(inNCLon)) == 3) {
    x <- as.vector(inNCLon[, ncol(inNCLon):1,])
    y <- as.vector(inNCLat[, ncol(inNCLat):1,])
  }

now it works. Because there was a similar problem in another program, I thought its the same.

[Schuch666]

Hi @bakamotokatas

No problem, there is no time wasted, the change from rgdal to sf was a pending task and since it is not possible to test all possibilities for the wrf_raster function all the input is welcome.

This last version of the wrf_raster contain some changes including your suggestion, an option to project the output to latlon, use_sf and a fix in layer name for 2d variables.

Let us know if you find ant new issues