Always necessary!

Let us say that you have a namelist.wps like this:

&share
 wrf_core = 'ARW',
 max_dom = 2,
 start_date = '2014-10-03_00:00:00','2014-10-03_00:00:00'
 end_date   = '2014-10-11_00:00:00','2014-10-11_00:00:00'
 interval_seconds = 21600
 io_form_geogrid = 2,
/

&geogrid
 parent_id         =   1,   1,
 parent_grid_ratio =   1,   3,
 i_parent_start    =   1,  35,
 j_parent_start    =   1,  33,
 e_we              =  95,  67,
 e_sn              =  90,  64,

 geog_data_res = 'usgs','usgs',
 dx = 9000,
 dy = 9000,
 map_proj = 'lambert',
 ref_lat   = -23.5,
 ref_lon   = -46.353519,
 truelat1  = -23.5,
 truelat2  = -23.5,
 stand_lon = -46.353519
 geog_data_path = '/PATH/WPS_GEOG'
/

&ungrib
 out_format = 'WPS',
 prefix = 'FILE',
/

&metgrid
 fg_name = 'FILE'
 io_form_metgrid = 2, 
/

Then you just do

./geogrid.exe

This will generate the files:

  • geo_em.d01.nc
  • geo_em.d02.nc

now let us plot the domains if they are in the right place!

First let us read the grid from the geo_em.d02.nc file, let us add the coatlines and also another shapefile for reference

library(eixport)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1
# d02
g <- wrf_grid("/media/sergio/ext51/WRF4/WPS/geo_em.d02.nc", type = "geo")
#> using grid info from: /media/sergio/ext51/WRF4/WPS/geo_em.d02.nc 
#> Number of lat points 125
#> Number of lon points 160
rmsp <- st_read("/media/sergio/0A9AD66165F337622/MEGA/PhDgrive/Dados/OD2007-METRO/od.shp")
#> Reading layer `od' from data source 
#>   `/media/sergio/0A9AD66165F337622/MEGA/PhDgrive/Dados/OD2007-METRO/od.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 460 features and 5 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -47.2086 ymin: -24.0642 xmax: -45.69495 ymax: -23.18363
#> Geodetic CRS:  WGS 84

plot(g$geometry,axes = T, reset = F, lwd = 0.1)
plot(rmsp$geometry, add= T)

And then with leaflet would be

library(leaflet)

#mapview(g)
leaflet::leaflet(g$geometry) %>%
  leaflet::addPolylines(color = "black", weight = 0.5) %>%
  leaflet::addTiles()