2

I try to convert Lambert conformal coordinates to lat/lon (WGS84) and I have used wgrib2, but the result is biased.

Command:

wgrib2 "mypath" -match "10m...." -new_grid_winds grid -new_grid_interpolation neighbor -new_grid latlon 108:129:0.25 16:65:0.25 "outputpath"

results with:

enter image description here

while it should be like that (from windy.com)

enter image description here

grib file:

Grib2 file

Grib2json file

Frank Liao
  • 855
  • 1
  • 8
  • 25
  • Can you share the grib file? – msi_gerva Sep 13 '18 at 08:04
  • https://drive.google.com/file/d/1eJL4TEnnVHkrLSaXrFPJxTYhmOyoN1Ig/view?usp=sharing – Frank Liao Sep 13 '18 at 08:19
  • Can you give exact command how you interpolate to the new values using wgrib2? I am not familiar with it, but I tried just to convert to the netCDF using `wgrib2 infile -netcdf outfile` and checked the result by plotting wind vectors with Python and the results looked good to me. – msi_gerva Sep 13 '18 at 13:53
  • I use wgrib2 convert to latlon, then I use grib2json convert it to json file and then I use /leaflet-velocity to show https://github.com/danwild/leaflet-velocity – Frank Liao Sep 14 '18 at 02:17
  • I was hoping for an exact command how you converted the file. Can you share the command and also the output from grib2json? As I said, if I converted to netCDF, made a plot, everything was ok... – msi_gerva Sep 14 '18 at 08:42
  • grib2json --names --data --output output.json mygrib.grib – Frank Liao Sep 14 '18 at 09:10
  • https://drive.google.com/file/d/1Pel5rwZ4HXFCF3UkpvCkhrVD-86Ko149/view?usp=sharing – Frank Liao Sep 14 '18 at 09:10
  • I was thinking about the exact command in the question. Is it really like this: `wgrib2 "mypath" -match "10m...." -new_grid_winds grid -new_grid_interpolation neighbor -new_grid latlon 108:129:0.25 16:65:0.25 "outputpath"`? If yes, then it did not work for me as I got messages: `**** ERROR: local table = 0 is not allowed, set to 1 ***`. I used command like this:`wgrib2 M-A0064-000.grb2 -match "10m" -new_grid_winds grid -new_grid_interpolation neighbor -new_grid latlon 108:129:0.25 16:65:0.25 newtest.nc` – msi_gerva Sep 14 '18 at 09:40
  • I got the same message... "**** ERROR: local table = 0 is not allowed, set to 1 ***" – Frank Liao Sep 17 '18 at 00:48

2 Answers2

2

I think there might be some flaws in the initial grib file. I converted the grib file to netCDF using wgrib2 and after that made some plots using Python and the result is not good.

Thing is, when I make the plot of temperature and overlay that with wind vectors, it looks ok. Problem is, when I also add the coastline, I see that the location of Taiwan island and also the main continent does not match with coastline drawn from the database.

Therefore I assume there is something bad in the initial gribfile - either the coordinates (start and endpoint or the step) are not very good and the coordinates written to the netCDF are not correct.

My code is here, if interested:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
import json
# -------------------------------
# read the json file:
with open('2018091312.json','r') as f:
    data = json.load(f)
# -------------------------------
lo1,lo2,la1,la2 = 108,142.75,16,23.75
dx,dy=0.25,0.25
nx,ny=140,32
udata=np.array(data[0]['data'],dtype='float32');udata=np.reshape(udata,(ny,nx));
vdata=np.array(data[1]['data'],dtype='float32');vdata=np.reshape(vdata,(ny,nx));
londata=np.arange(lo1,lo2+dx,dx);
latdata=np.arange(la1,la2+dy,dy);
londata,latdata=np.meshgrid(londata,latdata)
# -------------------------------
# -------------------------------
ncin=Dataset('test.nc');
lons=ncin.variables['longitude'][:];
lats=ncin.variables['latitude'][:];
u10=np.squeeze(ncin.variables['UGRD_10maboveground'][:])
v10=np.squeeze(ncin.variables['VGRD_10maboveground'][:])
t2=np.squeeze(ncin.variables['TMP_surface'][:])
ncin.close();
# -------------------------------
xlim=(np.min(lons),np.max(lons));
ylim=(np.min(lats),np.max(lats));
# -------------------------------
plt.figure(figsize=(8, 8))
m = Basemap(projection='cyl', resolution='i',
            llcrnrlat=ylim[0], urcrnrlat=ylim[1],
            llcrnrlon=xlim[0], urcrnrlon=xlim[1], )
xx,yy=m(lons,lats);
m.pcolormesh(lons,lats,t2,vmin=273.,vmax=300.);
skipx=skipy=16
m.quiver(xx[::skipy,::skipx],yy[::skipy,::skipx],u10[::skipy,::skipx],v10[::skipy,::skipx],scale=20.0,units='inches');
# ------------------------------------------
plt.savefig('test_withoutland.png',bbox_inches='tight')
m.drawcoastlines()
m.drawlsmask(land_color = "#ddaa66")
plt.savefig('test_withland.png',bbox_inches='tight')
plt.show()
# ------------------------------------------
skipx,skipy=2,2
plt.figure(figsize=(8, 8))
m = Basemap(projection='cyl', resolution='i',
            llcrnrlat=ylim[0], urcrnrlat=ylim[1],
            llcrnrlon=xlim[0], urcrnrlon=xlim[1], )
xx,yy=m(londata,latdata);
m.pcolormesh(lons,lats,t2,vmin=273.,vmax=300.);
m.quiver(xx[::skipy,::skipx],yy[::skipy,::skipx],udata[::skipy,::skipx],vdata[::skipy,::skipx],scale=20.0,units='inches');
# ------------------------------------------
m.drawcoastlines()
m.drawlsmask(land_color = "#ddaa66")
plt.savefig('test_json.png',bbox_inches='tight')
plt.show()

And the result looks like this (the test with JSON file): enter image description here

The conversion from grib to newCDF, I did like this:

wgrib2 M-A0064-000.grb2 -netcdf test.nc
msi_gerva
  • 2,021
  • 3
  • 22
  • 28
0

There are some weird definitions in the WRF LCC that you need to keep in mind when doing your reprojections. This website (unaffiliated) details most of it using python.

https://fabienmaussion.info/2018/01/06/wrf-projection/

Heymans
  • 128
  • 1
  • 9