source /work/ROMO/anaconda3/bin/activate geo

DOM=${1?Provide a domain eg, 12US1 or HEMIS}

if [ ${DOM} == 12US1 ]; then
  PROJ4="+proj=lcc +lat_1=33.0 +lat_2=45.0 +lat_0=40.0 +lon_0=-97.0 +y_0=1728000.0 +x_0=2556000.0 +a=6370000.0 +b=6370000.0 +to_meter=12000.0 +no_defs"
  NX=459
  NY=299
  CLIP="-140 10 -40 60"
elif [ ${DOM} == 12US2 ]; then
  PROJ4="+proj=lcc +lat_1=33.0 +lat_2=45.0 +lat_0=40.0 +lon_0=-97.0 +y_0=1620000.0 +x_0=2412000.0 +a=6370000.0 +b=6370000.0 +to_meter=12000.0 +no_defs"
  NX=396
  NY=246
  CLIP="-140 10 -40 60"
elif [ ${DOM} == 36US3 ]; then
  PROJ4="+proj=lcc +lat_0=40 +lon_0=-97 +lat_1=33 +lat_2=45 +x_0=2952000 +y_0=2772000 +R=6370000 +to_meter=36000 +no_defs"
  NX=172
  NY=148
  CLIP="-150 5 -30 70"
elif [ ${DOM} == HEMIS || ${DOM} == 108NHEMI2 ]; then
  PROJ4="+proj=stere +lat_0=90.0 +lat_ts=45.0 +lon_0=-98.0 +y_0=10098000.0 +x_0=10098000.0 +a=6370000.0 +b=6370000.0 +to_meter=108000.0 +no_defs"
  NX=187
  NY=187
  CLIP="-180 -20 180 91"
  OPTS=" -wrapdateline ${OPTS}"
elif [ ${DOM} == 1188NHEMI2 ]; then
  PROJ4="+proj=stere +lat_0=90.0 +lat_ts=45.0 +lon_0=-98.0 +y_0=10098000.0 +x_0=10098000.0 +a=6370000.0 +b=6370000.0 +to_meter=1188000.0 +no_defs"
  NX=17
  NY=17
  CLIP="-180 -20 180 91"
  OPTS=" -wrapdateline ${OPTS}"
fi

ogr2ogr -skipfailures -dim XY -clipsrc ${CLIP} ${OPTS} -t_srs "${PROJ4}"  intermed/tz_world_hours_${DOM}.geojson intermed/tz_world_hours.geojson
gdal_rasterize -of netCDF -at -a hours  -ts ${NX} ${NY} -te 0 0 ${NX} ${NY} -l tz_world_hours intermed/tz_world_hours_${DOM}.geojson intermed/tz_world_hours_${DOM}.nc
# rm intermed/tz_world_hours_${DOM}.geojson