Nishadh KA

Working with WRF CHEM part4 cbe domain run

2014-07-13


  1. Coimbatore domain made with four nests has been run with WRF-EMS with resolution starting from 100>27>9>3>1km. The domain is made using dwiz application in WRF-EMS.
  2. However, to test the domain above the Coimbatore region, the output was subject for wrfncxnj.py utility. It took 4 hours and 33 minutes to run this domain in IBM X4 series server computer with 16 GB ram.
  3. To overlay the Coimbatore city shape file above the WRF output NetCDF from wrfcnxnj.py. The city shapefile of Coimbatore converted into custom projection as per this, this and this. The final code look likes this. But unfortunately the WRF output land far above north of city boundary so it has to edited to involve the Coimbatore region.
  4. The server IBM started and made a copy of the domain cbeh1 into cbeh2 and open this domain and edited the domain in dwiz to make the domain cover the Coimbatore region.
  5. To ascertain the domain coverage the geogrid formed by wrf ems can be visualized in QGIS platform, for this the NC file was subject to open in qgis, but it ends with an error saying not supported format.
  6. So the option is to convert this NetCDF file into geotiff and visualize it in the QGIS.
  7. One advantage observed was the shapefile converted into the defined projection is giving the area derived information in a square kilometer, this was problem in road network assessment in the Coimbatore region, so the road network has to be converted into this projection and problem can rectified.
  8. To convert NC file into geotiff, name of the variables were chosen by ncview. However, this ends up in error
  9. So made a last attempt of viewing nc file and its various variables from pythonic way, the command used was

    from osgeo import gdal
    from osgeo import osr
    import numpy
    import numpy.ma as ma
    datafile = 'geo_em.d04.nc'
    ds_in = gdal.Open(datafile)
    subdatasets = ds_in.GetSubDatasets()
    subdatasets
    ```
    these command gives the list of variables in the geo_em.d04.nc. In which the variables CLONG and CLAT was considered as lon and lat and variabnle GREENFRAC was used as parameter for geotiff conversion script. The script used were 
    ```
    from osgeo import gdal
    from osgeo import osr
    import numpy
    import numpy.ma as ma
    datafile = 'geo_em.d04.nc'
    proj_out = osr.SpatialReference()
    proj_out.SetMercator(0.0, 115.02, 0.98931892612652, 0.0, 0.0)
    ds_in = gdal.Open(datafile)
    #subdatasets = ds_in.GetSubDatasets()
    #variables = []
    #for subdataset in subdatasets:
    #     variables.append(subdataset[1].split(" ")[1])
    ds_lon = gdal.Open('NETCDF:"geo_em.d04.nc":CLONG')
    ds_lat = gdal.Open('NETCDF:"geo_em.d04.nc":CLAT')
    longs = ds_lon.GetRasterBand(1).ReadAsArray()
    lats = ds_lat.GetRasterBand(1).ReadAsArray()
    ds_lon = None
    ds_lat = None
    proj_gcp = osr.SpatialReference()
    proj_gcp.ImportFromEPSG(4326)
    transf = osr.CoordinateTransformation(proj_gcp, proj_out)
    ul = transf.TransformPoint(float(longs[0][0]), float(lats[0][0]))
    lr = transf.TransformPoint(float(longs[len(longs)-1][len(longs[0])-1]),
    float(lats[len(longs)-1][len(longs[0])-1]))
    ur = transf.TransformPoint(float(longs[0][len(longs[0])-1]),
    float(lats[0][len(longs[0])-1]))
    ll = transf.TransformPoint(float(longs[len(longs)-1][0]),
    float(lats[len(longs)-1][0]))
    gt0 = ul[0]
    gt1 = (ur[0] - gt0) / len(longs[0])
    gt2 = (lr[0] - gt0 - gt1*len(longs[0])) / len(longs)
    gt3 = ul[1]
    gt5 = (ll[1] - gt3) / len(longs)
    gt4 = (lr[1] - gt3 - gt5*len(longs) ) / len(longs[0])
    gt = (gt0,gt1,gt2,gt3,gt4,gt5)
    ds_u10 = gdal.Open('NETCDF:"'+datafile+'":GREENFRAC')
    #i=23
    #num_bands = ds_u10.RasterCount
    x_size = ds_u10.RasterXSize
    y_size = ds_u10.RasterYSize
    ds_p_b = ds_u10.GetRasterBand(1).ReadAsArray()
    driver = gdal.GetDriverByName( 'GTiff' )
    ds_out = driver.Create( 'geo_em.d04.tiff', x_size, y_size, 1, gdal.GDT_Float32 )
    ds_out.SetGeoTransform( gt )
    ds_out.SetProjection(proj_out.ExportToWkt())
    ds_out.GetRasterBand(1).WriteArray( ds_p_b )
    

    By this script the geotiff file geo_em.d04.tiff was created and it gets cover Coimbatore city boundary. The wrf -ems domain made sure that it involves the Coimbatore region extent.

  10. Made a copy of the cbe_h1 folder in runs folder of WRF ems and kept it to use for input.wps and input.real in wrf chem simulation over Coimbatore domain