changliao1025/exam_culled_mpas_mesh_cell.py

## exam_culled_mpas_mesh_cell.py
#this function can be used to extract the cell information from the MPAS mesh file and check whether the mesh cell is valid or not
#Author: Chang Liao, chang.liao@pnnl.gov
import os, sys
import math
import importlib.util
import numpy as np
import netCDF4 as nc
from osgeo import gdal, osr, ogr
iFlag_cython = importlib.util.find_spec("cython")
#convert from mpas 360 to 180 degree for longitude
if iFlag_cython is not None:
    from pyflowline.algorithms.cython.kernel import convert_360_to_180
else:
    from pyearth.gis.geometry.convert_longitude_range import convert_360_to_180
gdal.UseExceptions()

def plot_wkt_polygon(pPolygon_wkt, sFilename_png_out):
    #this is optional, if you want to plot the polygon, you can use this function, it calls some functions from the pyearth package
    import matplotlib.pyplot as plt
    from matplotlib.patches import Polygon
    from matplotlib.collections import PatchCollection
    import matplotlib as mpl
    import cartopy.crs as ccrs
    from cartopy.io.img_tiles import OSM
    from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
    from pyearth.gis.location.get_geometry_coordinates import get_geometry_coordinates
    from pyearth.visual.map.map_servers import calculate_zoom_level, calculate_scale_denominator
    from pyearth.visual.map.zebra_frame import zebra_frame
    pSRS_wgs84 = ccrs.PlateCarree()  # for latlon data only
    #get the list of coordinates in the wktpolygon
    pPolygon = ogr.CreateGeometryFromWkt(pPolygon_wkt)
    aCoords_gcs = get_geometry_coordinates(pPolygon)
    aLon = aCoords_gcs[:,0]
    aLat = aCoords_gcs[:,1]
    dLon_max = np.max(aLon)
    dLon_min = np.min(aLon)
    dLat_max = np.max(aLat)
    dLat_min = np.min(aLat)
    pProjection_map = ccrs.Orthographic(central_longitude=0.50*(
            dLon_max+dLon_min),  central_latitude=0.50*(dLat_max+dLat_min), globe=None)

    fig = plt.figure(dpi=150)
    fig.set_figwidth(8)
    fig.set_figheight(8)
    ax = fig.add_axes([0.08, 0.1, 0.62, 0.7], projection=pProjection_map)
    image_size = [1000, 1000]
    aExtent = [dLon_min, dLon_max, dLat_min, dLat_max]
    scale_denominator = calculate_scale_denominator(aExtent, image_size)
    pSrc = osr.SpatialReference()
    pSrc.ImportFromEPSG(3857) # mercator
    pProjection = pSrc.ExportToWkt()
    iBasemap_zoom_level = calculate_zoom_level(scale_denominator, pProjection)
    osm_tiles = OSM()
    #Add the OSM image to the map
    ax.add_image(osm_tiles, iBasemap_zoom_level)
    aPolygon = list()
    aPolygon.append(aCoords_gcs)
    aPatch = [Polygon(poly, closed=True) for poly in aPolygon]
    pPC = PatchCollection(aPatch, alpha=0.5,
                                  edgecolor='black',
                                  facecolor='none',
                                  linewidths=2,
                                  transform=pSRS_wgs84)
    ax.add_collection(pPC)
    marginx = (dLon_max - dLon_min) / 20
    marginy = (dLat_max - dLat_min) / 20
    if (dLat_max + marginy)> 90:
        dLat_max = 90
    else:
        dLat_max = dLat_max + marginy
    if (dLat_min - marginy) < -90:
        dLat_min = -90
    else:
        dLat_min = dLat_min - marginy
    if (dLon_max + marginx) > 180:
        dLon_max = 180
    else:
        dLon_max = dLon_max + marginx
    if (dLon_min - marginx) < -180:
        dLon_min = -180
    else:
        dLon_min = dLon_min - marginx
    aExtent_extend = [dLon_min, dLon_max, dLat_min, dLat_max]
    minx, maxx, miny, maxy = aExtent_extend
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
                      linewidth=1, color='gray', alpha=0.5, linestyle='--',
                      xlocs=np.arange(minx, maxx+(maxx-minx)/9, (maxx-minx)/8),
                      ylocs=np.arange(miny, maxy+(maxy-miny)/9, (maxy-miny)/8))
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlocator = mpl.ticker.MaxNLocator(4)
    gl.ylocator = mpl.ticker.MaxNLocator(4)
    gl.xlabel_style = {'size': 10, 'color': 'k', 'rotation': 0, 'ha': 'right'}
    gl.ylabel_style = {'size': 10, 'color': 'k',
                       'rotation': 90, 'weight': 'normal'}

    ax.set_extent(aExtent_extend, crs = pSRS_wgs84)
    ax.zebra_frame(crs=pSRS_wgs84, iFlag_outer_frame_in=1)
    plt.savefig(sFilename_png_out, bbox_inches='tight')
    print('The plot is saved to: ', sFilename_png_out)
    return

def check_mesh_cell_polygon(sFilename_mpas_mesh_netcdf_culled, aCellID_in,
                             iFlag_plot_in = None ,
                            sFilename_mpas_mesh_netcdf_base = None):
    #sFilename_mpas_mesh_netcdf_culled, the culled mesh cell netcdf
    #aCellID_in, the list of cell ID to be checked
    #iFlag_plot_in, if it is not None, it will plot the cell polygon
    #sFilename_mpas_mesh_netcdf_base, the base mesh cell netcdf, it will check whether the base cell is valid or not, using mesh cell center
    #first, we need to read the mesh file and copy it to the new file
    pDatasets_mesh = nc.Dataset(sFilename_mpas_mesh_netcdf_culled, 'r')
    #get netcdf format
    format = pDatasets_mesh.file_format
    print('Netcdf format:', format)
    #read new netcdf
    for sKey, aValue in pDatasets_mesh.variables.items():
        if sKey == 'lonCell':
            lonCell0 = aValue

        if sKey == 'latCell':
            latCell0 = aValue

        if sKey == 'verticesOnCell':
            verticesOnCell0 = aValue

        if sKey == 'indexToCellID':
            indexToCellID0 = aValue

        if sKey == 'lonVertex':
            lonVertex0 = aValue

        if sKey == 'latVertex':
            latVertex0 = aValue

    aLongitudeCell = lonCell0[:] / math.pi * 180
    aLatitudeCell = latCell0[:] / math.pi * 180
    aLongitudeVertex = lonVertex0[:] / math.pi * 180
    aLatitudeVertex = latVertex0[:] / math.pi * 180
    aVertexOnCell = verticesOnCell0[:]
    aIndexToCellID = indexToCellID0[:]
    ncell = len(aIndexToCellID)
    aIndexToCellID = np.array(aIndexToCellID)

    if sFilename_mpas_mesh_netcdf_base is not None:
        pDatasets_mesh_base = nc.Dataset(sFilename_mpas_mesh_netcdf_base, 'r')
        for sKey, aValue in pDatasets_mesh_base.variables.items():
            if sKey == 'lonCell':
                lonCell1 = aValue

            if sKey == 'latCell':
                latCell1 = aValue

            if sKey == 'verticesOnCell':
                verticesOnCell1 = aValue

            if sKey == 'indexToCellID':
                indexToCellID1 = aValue

            if sKey == 'lonVertex':
                lonVertex1 = aValue

            if sKey == 'latVertex':
                latVertex1 = aValue

        aLongitudeCell_base = lonCell1[:] / math.pi * 180
        aLatitudeCell_base = latCell1[:] / math.pi * 180
        aLongitudeCell_base = np.array(aLongitudeCell_base)
        aLatitudeCell_base = np.array(aLatitudeCell_base)
        aLongitudeVertex_base = lonVertex1[:] / math.pi * 180
        aLatitudeVertex_base = latVertex1[:] / math.pi * 180
        aVertexOnCell_base = verticesOnCell1[:]
        aIndexToCellID_base = indexToCellID1[:]
        ncell_base = len(aIndexToCellID_base)
        aIndexToCellID_base = np.array(aIndexToCellID_base)


    for lCellID_in in aCellID_in:
        if lCellID_in > ncell or lCellID_in < 0:
            print('Error: the cell ID is out of range')
            return

        #use the cell ID to obtain its index, all MPAS index starts from 1!

        lCellIndex = np.where(aIndexToCellID == lCellID_in)
        dLon_center_360 = float(aLongitudeCell[lCellIndex])
        dLon_center = convert_360_to_180 (dLon_center_360)
        dLat_center =  float(aLatitudeCell[lCellIndex])
        aVertexOnCellIndex = np.array(aVertexOnCell[lCellIndex,:])
        dummy0 = np.where(aVertexOnCellIndex > 0)
        aVertexIndex = aVertexOnCellIndex[dummy0] - 1
        aLonVertex = aLongitudeVertex[aVertexIndex]
        aLatVertex = aLatitudeVertex[aVertexIndex]
        nVertex = len(aLonVertex)
        if nVertex < 3:
            print("Vertex number is: ", nVertex)
            return
        else:
            ring = ogr.Geometry(ogr.wkbLinearRing)
            aCoords_gcs = np.full((nVertex,2), -9999.0, dtype=float)
            for j in range(nVertex):
                x1 = convert_360_to_180(aLonVertex[j])
                y1 = float(aLatVertex[j])
                ring.AddPoint(x1, y1)
                aCoords_gcs[j,0] = x1
                aCoords_gcs[j,1] = y1
                pass

            #add the closing point, which is the first vertex
            x1 = convert_360_to_180(aLonVertex[0])
            y1 = aLatVertex[0]
            ring.AddPoint(x1, y1) #double check
            pPolygon = ogr.Geometry(ogr.wkbPolygon)
            pPolygon.AddGeometry(ring)
            # Validate the geometry
            pPolygon_wkt = pPolygon.ExportToWkt()
            if not pPolygon.IsValid():
                print("Polygon is invalid...")
                #pPolygon = pPolygon.Buffer(0)  # Buffering by 0 can fix some geometry issues
                #pPolygon = pPolygon.Simplify(0.0001)  # Simplify with a tolerance
            else:
                print("Polygon is valid...")

            print("Polygon WKT:", pPolygon_wkt)
            if iFlag_plot_in is not None:
                #use a simple method to plot the polygon
                #get the path where this script is located
                sPath_script = os.path.dirname(os.path.abspath(__file__))
                sFilename_png = 'mesh_cell_' + str(lCellID_in) + '_culled.png'
                sFilename_png_out = os.path.join(sPath_script, sFilename_png)
                plot_wkt_polygon(pPolygon_wkt, sFilename_png_out)
                if sFilename_mpas_mesh_netcdf_base is not None:
                    #find the base mesh cell using the center lon and lat
                    lCell_index_base = np.where((aLongitudeCell_base == dLon_center_360) & (aLatitudeCell_base == dLat_center))
                    #check we can find a base cell or not
                    if len(lCell_index_base) == 1 and lCell_index_base[0] != 0:
                        lCellID_base = aIndexToCellID_base[lCell_index_base]
                        aVertexOnCellIndex = np.array(aVertexOnCell_base[lCell_index_base,:])
                        dummy0 = np.where(aVertexOnCellIndex > 0)
                        aVertexIndex = aVertexOnCellIndex[dummy0] - 1
                        aLonVertex_base = aLongitudeVertex_base[aVertexIndex]
                        aLatVertex_base = aLatitudeVertex_base[aVertexIndex]
                        nVertex = len(aLonVertex)
                        if nVertex < 3:
                            print("Base vertex number is: ", nVertex)
                            return
                        else:
                            ring = ogr.Geometry(ogr.wkbLinearRing)
                            aCoords_gcs = np.full((nVertex,2), -9999.0, dtype=float)
                            for j in range(nVertex):
                                x1 = convert_360_to_180(aLonVertex_base[j])
                                y1 = float(aLatVertex_base[j])
                                ring.AddPoint(x1, y1)
                                aCoords_gcs[j,0] = x1
                                aCoords_gcs[j,1] = y1
                                pass

                            #add the closing point, which is the first vertex
                            x1 = convert_360_to_180(aLonVertex_base[0])
                            y1 = aLatVertex_base[0]
                            ring.AddPoint(x1, y1) #double check
                            pPolygon = ogr.Geometry(ogr.wkbPolygon)
                            pPolygon.AddGeometry(ring)
                            # Validate the geometry
                            pPolygon_wkt_base = pPolygon.ExportToWkt()
                            if not pPolygon.IsValid():
                                print("Base polygon is invalid...")
                                #pPolygon = pPolygon.Buffer(0)  # Buffering by 0 can fix some geometry issues
                                #pPolygon = pPolygon.Simplify(0.0001)  # Simplify with a tolerance
                            else:
                                print("Base polygon is valid...")

                            sFilename_png = 'mesh_cell_' + str(lCellID_in) + '_base.png'
                            sFilename_png_out = os.path.join(sPath_script, sFilename_png)
                            plot_wkt_polygon(pPolygon_wkt_base, sFilename_png_out)
                    else:
                        print("The base cell is not found...")
                        continue
                pass


    return

#create a main call
if __name__ == '__main__':
    sFilename_mpas_mesh_netcdf_culled = '/compyfs/liao313/04model/pyflowline/conus/pyflowline20241201019/jigsaw/out/invert_mesh.nc'
    sFilename_mpas_mesh_netcdf_base = '/compyfs/liao313/04model/pyflowline/conus/pyflowline20241201019/jigsaw/out/base_mesh.nc'
    aCellID_in = list()
    #Polygon is invalid... 1672
    #Polygon WKT: POLYGON ((122.280664576308 18.2221473603637 0,122.262233339828 18.2160835887451 0,122.263243646392 18.1964473296003 0,122.279796144809 18.1912927994405 0,122.287660730398 18.2094542663177 0,122.278850506809 18.2232593941033 0,122.280664576308 18.2221473603637 0))
    #Polygon is invalid... 13994
    #Polygon WKT: POLYGON ((-63.9008296243521 50.3336863925212 0,-63.9273355026816 50.3213375919327 0,-63.9062274615777 50.3021870895725 0,-63.9271251813239 50.3195213928764 0,-63.8814145262508 50.3029476266043 0,-63.8690894249787 50.311717913351 0,-63.8704197044186 50.3233279141183 0,-63.9008296243521 50.3336863925212 0))
    #Polygon is invalid... 15428
    #Polygon WKT: POLYGON ((-48.7470719201964 61.4094980978841 0,-48.7078667093716 61.4115580717041 0,-48.7076968290549 61.4159935394359 0,-48.7244691494012 61.4238440398355 0,-48.7482711956022 61.4177357340897 0,-48.7500637709318 61.4113714878841 0,-48.7318223408571 61.407312966789 0,-48.7470719201964 61.4094980978841 0))
    #Polygon is invalid... 17062
    #Polygon WKT: POLYGON ((-129.539599690119 70.1323831312393 0,-129.535795196565 70.124261028838 0,-129.512483903554 70.1172844555891 0,-129.465730846556 70.1172188707654 0,-129.459111064388 70.1364237439907 0,-129.465652367021 70.1172520619199 0,-129.498714599157 70.149964048113 0,-129.539599690119 70.1323831312393 0))
    #Polygon is invalid... 17793
    #Polygon WKT: POLYGON ((-71.4446320396707 10.4020008778489 0,-71.4516631039428 10.3961118391854 0,-71.456830793219 10.4214327874876 0,-71.4760303796148 10.4223603510441 0,-71.4842067505045 10.4056232943428 0,-71.4713271460859 10.3935027872622 0,-71.4446320396707 10.4020008778489 0))
    #Polygon is invalid... 24092
    #Polygon WKT: POLYGON ((-15.5382028906378 11.5876827405335 0,-15.535595747457 11.5850973512438 0,-15.5099707980737 11.591618233233 0,-15.5071308491261 11.599961078997 0,-15.5246821263325 11.613672338718 0,-15.5098012098433 11.6054745495611 0,-15.5420411080539 11.6092003139624 0,-15.5382028906378 11.5876827405335 0))
    #Polygon is invalid... 24513
    #Polygon WKT: POLYGON ((-43.3410601521804 -2.35730294671818 0,-43.3101695695966 -2.36488149251201 0,-43.3081346848379 -2.35641950002847 0,-43.3166767565955 -2.33487861443757 0,-43.3342404233832 -2.3423723030464 0,-43.3405705233804 -2.36152207379487 0,-43.3410601521804 -2.35730294671818 0))
    #Polygon is invalid... 33144
    #Polygon WKT: POLYGON ((119.364514438706 -1.17964356953732 0,119.348411574817 -1.17792521472226 0,119.339621749198 -1.18809595090908 0,119.34848087937 -1.20993515223943 0,119.34738302738 -1.2089944663736 0,119.365480682366 -1.20127151411875 0,119.364514438706 -1.17964356953732 0))
    #Polygon is invalid... 35167
    #Polygon WKT: POLYGON ((-35.9944101244893 -5.05398122354815 0,-35.9801257279569 -5.06874982595671 0,-35.9656043467338 -5.06703324195992 0,-35.9708716349697 -5.04776250575661 0,-35.9648702774745 -5.06163712690284 0,-35.9884656481591 -5.03975398282233 0,-35.9944101244893 -5.05398122354815 0))
    #Polygon is invalid... 35616
    #Polygon WKT: POLYGON ((-157.006243393183 71.2109730489155 0,-157.010511724013 71.1907258531127 0,-156.984279984155 71.1778101013443 0,-156.986514943451 71.1781925881701 0,-156.93530954857 71.1846969871186 0,-156.948468763097 71.2107407300836 0,-157.006243393183 71.2109730489155 0))
    #Polygon is invalid... 35641
    #Polygon WKT: POLYGON ((33.7201130735598 28.0767992360369 0,33.7389806225413 28.0809644371193 0,33.734809560818 28.0775273580535 0,33.734179155333 28.0972618198086 0,33.7132054411147 28.1016558473581 0,33.7023215895678 28.086756860363 0,33.7201130735598 28.0767992360369 0))
    #Polygon is invalid... 38117
    #Polygon WKT: POLYGON ((148.343510865271 -5.59033620753651 0,148.355315036408 -5.57459782311827 0,148.3452916253 -5.55419706603241 0,148.343320630905 -5.55389792166596 0,148.329597677703 -5.5616335416741 0,148.326235132406 -5.57660836257517 0,148.345384562197 -5.59062056647374 0,148.343510865271 -5.59033620753651 0))
    #Polygon is invalid... 41549
    #Polygon WKT: POLYGON ((137.74958112273 -35.1129190969386 0,137.728498752014 -35.1157816659266 0,137.763899893261 -35.0954158510455 0,137.735133286168 -35.0906702014205 0,137.719432127827 -35.1031570555177 0,137.722144962911 -35.1135022577956 0,137.74958112273 -35.1129190969386 0))
    #Polygon is invalid... 41660
    #Polygon WKT: POLYGON ((-16.1689593883591 12.2942482599369 0,-16.1656619757907 12.2913379381728 0,-16.1501863914894 12.2909515078012 0,-16.1388797897555 12.302421020704 0,-16.1457460487068 12.3239541811496 0,-16.1435360602088 12.3216513216158 0,-16.1733726860637 12.3104519598962 0,-16.1689593883591 12.2942482599369 0))
    #Polygon is invalid... 42129
    #Polygon WKT: POLYGON ((-71.873745052447 10.2041331098005 0,-71.8695608031783 10.2010706472204 0,-71.8505573208394 10.2005214388756 0,-71.8431490329504 10.2148235313675 0,-71.8501723767827 10.2314272143016 0,-71.8495997694914 10.2310373849819 0,-71.8720248298625 10.2204736226489 0,-71.873745052447 10.2041331098005 0))
    #Polygon is invalid... 45022
    #Polygon WKT: POLYGON ((-3.49990419525511 50.4458670833471 0,-3.50109363562478 50.4650696303941 0,-3.53783879506921 50.4773819921189 0,-3.52408079611428 50.4765892264093 0,-3.54618447523524 50.4521880445044 0,-3.52947657697973 50.4430455316903 0,-3.49990419525511 50.4458670833471 0))
    #Polygon is invalid... 46969
    #Polygon WKT: POLYGON ((-42.8889909455578 61.0765107756435 0,-42.9121083332119 61.0772667813284 0,-42.916511986311 61.0595397260809 0,-42.8723073974284 61.0505799864484 0,-42.869995504361 61.0608183805887 0,-42.8707395979006 61.0514417996377 0,-42.8889909455578 61.0765107756435 0))
    #Polygon is invalid... 50195
    #Polygon WKT: POLYGON ((-94.368666425969 68.9908495945327 0,-94.3175518357232 68.9953417755255 0,-94.3815668165992 68.9736632134828 0,-94.3425952667274 68.9661781569792 0,-94.2997990225276 68.9815048216357 0,-94.313890214263 68.9949962318141 0,-94.368666425969 68.9908495945327 0))
    #Polygon is invalid... 51427
    #Polygon WKT: POLYGON ((-91.9305657042177 29.8364298178733 0,-91.9471611860289 29.8502110560581 0,-91.9683012276131 29.8435993051381 0,-91.9705116516197 29.826359785042 0,-91.966002264764 29.8170637700613 0,-91.9592633280558 29.8132609164735 0,-91.931947199102 29.8332217173152 0,-91.9315475545814 29.8412767907539 0,-91.9305657042177 29.8364298178733 0))
    #Polygon is invalid... 52022
    #Polygon WKT: POLYGON ((119.452814084782 26.0049533277492 0,119.461850081725 25.9959711700541 0,119.484939468505 25.9981496514189 0,119.490311905252 26.0107563312908 0,119.475593022448 26.0293297088006 0,119.481284836318 26.0265138358557 0,119.452814084782 26.0049533277492 0))
    #Polygon is invalid... 53684
    #Polygon WKT: POLYGON ((134.880996223007 34.6996946747361 0,134.885760101184 34.7052291283607 0,134.86918556633 34.7273923580406 0,134.846823160703 34.7248603694141 0,134.846104671205 34.7037126681144 0,134.84636206586 34.724189250281 0,134.861751304972 34.6951355631687 0,134.880996223007 34.6996946747361 0))
    #Polygon is invalid... 53873
    #Polygon WKT: POLYGON ((67.8621346406756 71.4536132089809 0,67.8603493762603 71.4692951262031 0,67.820678200099 71.4835092147059 0,67.7810267313965 71.4725449644734 0,67.7937258331278 71.4477814842135 0,67.7825844661424 71.4529756746679 0,67.8476571359303 71.4481797847844 0,67.8621346406756 71.4536132089809 0))
    #Polygon is invalid... 54221
    #Polygon WKT: POLYGON ((-42.3203661613626 62.709446071854 0,-42.3313270148944 62.7199910936606 0,-42.3741515666778 62.7235736362693 0,-42.3760364767373 62.7051599709813 0,-42.3803038139625 62.7147640854633 0,-42.3462106371608 62.6948465209412 0,-42.3274484550065 62.6974134802222 0,-42.3203661613626 62.709446071854 0))
    #Polygon is invalid... 55261
    #Polygon WKT: POLYGON ((-96.2548737107388 68.5209704407071 0,-96.2977960367591 68.5145042550792 0,-96.3095326985719 68.4931769995066 0,-96.2957756974138 68.4880380626558 0,-96.2514869022441 68.4879607107365 0,-96.2306622412806 68.5028910535988 0,-96.2581634506687 68.5217217190116 0,-96.2548737107388 68.5209704407071 0))
    #Polygon is invalid... 57712
    #aCellID_in.append(1672+1)
    #aCellID_in.append(13994+1)
    #aCellID_in.append(15428+1)
    aCellID_in.append(17062+1)
    #aCellID_in.append(17793+1)
    #aCellID_in.append(24092+1)
    #aCellID_in.append(24513+1)
    #aCellID_in.append(33144+1)
    #aCellID_in.append(35167+1)

    check_mesh_cell_polygon(sFilename_mpas_mesh_netcdf_culled, aCellID_in,  iFlag_plot_in = 1, sFilename_mpas_mesh_netcdf_base = sFilename_mpas_mesh_netcdf_base )
	#this function can be used to extract the cell information from the MPAS mesh file and check whether the mesh cell is valid or not
	#Author: Chang Liao, chang.liao@pnnl.gov
	import os, sys
	import math
	import importlib.util
	import numpy as np
	import netCDF4 as nc
	from osgeo import gdal, osr, ogr
	iFlag_cython = importlib.util.find_spec("cython")
	#convert from mpas 360 to 180 degree for longitude
	if iFlag_cython is not None:
	from pyflowline.algorithms.cython.kernel import convert_360_to_180
	else:
	from pyearth.gis.geometry.convert_longitude_range import convert_360_to_180
	gdal.UseExceptions()

	def plot_wkt_polygon(pPolygon_wkt, sFilename_png_out):
	#this is optional, if you want to plot the polygon, you can use this function, it calls some functions from the pyearth package
	import matplotlib.pyplot as plt
	from matplotlib.patches import Polygon
	from matplotlib.collections import PatchCollection
	import matplotlib as mpl
	import cartopy.crs as ccrs
	from cartopy.io.img_tiles import OSM
	from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
	from pyearth.gis.location.get_geometry_coordinates import get_geometry_coordinates
	from pyearth.visual.map.map_servers import calculate_zoom_level, calculate_scale_denominator
	from pyearth.visual.map.zebra_frame import zebra_frame
	pSRS_wgs84 = ccrs.PlateCarree() # for latlon data only
	#get the list of coordinates in the wktpolygon
	pPolygon = ogr.CreateGeometryFromWkt(pPolygon_wkt)
	aCoords_gcs = get_geometry_coordinates(pPolygon)
	aLon = aCoords_gcs[:,0]
	aLat = aCoords_gcs[:,1]
	dLon_max = np.max(aLon)
	dLon_min = np.min(aLon)
	dLat_max = np.max(aLat)
	dLat_min = np.min(aLat)
	pProjection_map = ccrs.Orthographic(central_longitude=0.50*(
	dLon_max+dLon_min), central_latitude=0.50*(dLat_max+dLat_min), globe=None)

	fig = plt.figure(dpi=150)
	fig.set_figwidth(8)
	fig.set_figheight(8)
	ax = fig.add_axes([0.08, 0.1, 0.62, 0.7], projection=pProjection_map)
	image_size = [1000, 1000]
	aExtent = [dLon_min, dLon_max, dLat_min, dLat_max]
	scale_denominator = calculate_scale_denominator(aExtent, image_size)
	pSrc = osr.SpatialReference()
	pSrc.ImportFromEPSG(3857) # mercator
	pProjection = pSrc.ExportToWkt()
	iBasemap_zoom_level = calculate_zoom_level(scale_denominator, pProjection)
	osm_tiles = OSM()
	#Add the OSM image to the map
	ax.add_image(osm_tiles, iBasemap_zoom_level)
	aPolygon = list()
	aPolygon.append(aCoords_gcs)
	aPatch = [Polygon(poly, closed=True) for poly in aPolygon]
	pPC = PatchCollection(aPatch, alpha=0.5,
	edgecolor='black',
	facecolor='none',
	linewidths=2,
	transform=pSRS_wgs84)
	ax.add_collection(pPC)
	marginx = (dLon_max - dLon_min) / 20
	marginy = (dLat_max - dLat_min) / 20
	if (dLat_max + marginy)> 90:
	dLat_max = 90
	else:
	dLat_max = dLat_max + marginy
	if (dLat_min - marginy) < -90:
	dLat_min = -90
	else:
	dLat_min = dLat_min - marginy
	if (dLon_max + marginx) > 180:
	dLon_max = 180
	else:
	dLon_max = dLon_max + marginx
	if (dLon_min - marginx) < -180:
	dLon_min = -180
	else:
	dLon_min = dLon_min - marginx
	aExtent_extend = [dLon_min, dLon_max, dLat_min, dLat_max]
	minx, maxx, miny, maxy = aExtent_extend
	gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
	linewidth=1, color='gray', alpha=0.5, linestyle='--',
	xlocs=np.arange(minx, maxx+(maxx-minx)/9, (maxx-minx)/8),
	ylocs=np.arange(miny, maxy+(maxy-miny)/9, (maxy-miny)/8))
	gl.xformatter = LONGITUDE_FORMATTER
	gl.yformatter = LATITUDE_FORMATTER
	gl.xlocator = mpl.ticker.MaxNLocator(4)
	gl.ylocator = mpl.ticker.MaxNLocator(4)
	gl.xlabel_style = {'size': 10, 'color': 'k', 'rotation': 0, 'ha': 'right'}
	gl.ylabel_style = {'size': 10, 'color': 'k',
	'rotation': 90, 'weight': 'normal'}

	ax.set_extent(aExtent_extend, crs = pSRS_wgs84)
	ax.zebra_frame(crs=pSRS_wgs84, iFlag_outer_frame_in=1)
	plt.savefig(sFilename_png_out, bbox_inches='tight')
	print('The plot is saved to: ', sFilename_png_out)
	return

	def check_mesh_cell_polygon(sFilename_mpas_mesh_netcdf_culled, aCellID_in,
	iFlag_plot_in = None ,
	sFilename_mpas_mesh_netcdf_base = None):
	#sFilename_mpas_mesh_netcdf_culled, the culled mesh cell netcdf
	#aCellID_in, the list of cell ID to be checked
	#iFlag_plot_in, if it is not None, it will plot the cell polygon
	#sFilename_mpas_mesh_netcdf_base, the base mesh cell netcdf, it will check whether the base cell is valid or not, using mesh cell center
	#first, we need to read the mesh file and copy it to the new file
	pDatasets_mesh = nc.Dataset(sFilename_mpas_mesh_netcdf_culled, 'r')
	#get netcdf format
	format = pDatasets_mesh.file_format
	print('Netcdf format:', format)
	#read new netcdf
	for sKey, aValue in pDatasets_mesh.variables.items():
	if sKey == 'lonCell':
	lonCell0 = aValue

	if sKey == 'latCell':
	latCell0 = aValue

	if sKey == 'verticesOnCell':
	verticesOnCell0 = aValue

	if sKey == 'indexToCellID':
	indexToCellID0 = aValue

	if sKey == 'lonVertex':
	lonVertex0 = aValue

	if sKey == 'latVertex':
	latVertex0 = aValue

	aLongitudeCell = lonCell0[:] / math.pi * 180
	aLatitudeCell = latCell0[:] / math.pi * 180
	aLongitudeVertex = lonVertex0[:] / math.pi * 180
	aLatitudeVertex = latVertex0[:] / math.pi * 180
	aVertexOnCell = verticesOnCell0[:]
	aIndexToCellID = indexToCellID0[:]
	ncell = len(aIndexToCellID)
	aIndexToCellID = np.array(aIndexToCellID)

	if sFilename_mpas_mesh_netcdf_base is not None:
	pDatasets_mesh_base = nc.Dataset(sFilename_mpas_mesh_netcdf_base, 'r')
	for sKey, aValue in pDatasets_mesh_base.variables.items():
	if sKey == 'lonCell':
	lonCell1 = aValue

	if sKey == 'latCell':
	latCell1 = aValue

	if sKey == 'verticesOnCell':
	verticesOnCell1 = aValue

	if sKey == 'indexToCellID':
	indexToCellID1 = aValue

	if sKey == 'lonVertex':
	lonVertex1 = aValue

	if sKey == 'latVertex':
	latVertex1 = aValue

	aLongitudeCell_base = lonCell1[:] / math.pi * 180
	aLatitudeCell_base = latCell1[:] / math.pi * 180
	aLongitudeCell_base = np.array(aLongitudeCell_base)
	aLatitudeCell_base = np.array(aLatitudeCell_base)
	aLongitudeVertex_base = lonVertex1[:] / math.pi * 180
	aLatitudeVertex_base = latVertex1[:] / math.pi * 180
	aVertexOnCell_base = verticesOnCell1[:]
	aIndexToCellID_base = indexToCellID1[:]
	ncell_base = len(aIndexToCellID_base)
	aIndexToCellID_base = np.array(aIndexToCellID_base)


	for lCellID_in in aCellID_in:
	if lCellID_in > ncell or lCellID_in < 0:
	print('Error: the cell ID is out of range')
	return

	#use the cell ID to obtain its index, all MPAS index starts from 1!

	lCellIndex = np.where(aIndexToCellID == lCellID_in)
	dLon_center_360 = float(aLongitudeCell[lCellIndex])
	dLon_center = convert_360_to_180 (dLon_center_360)
	dLat_center = float(aLatitudeCell[lCellIndex])
	aVertexOnCellIndex = np.array(aVertexOnCell[lCellIndex,:])
	dummy0 = np.where(aVertexOnCellIndex > 0)
	aVertexIndex = aVertexOnCellIndex[dummy0] - 1
	aLonVertex = aLongitudeVertex[aVertexIndex]
	aLatVertex = aLatitudeVertex[aVertexIndex]
	nVertex = len(aLonVertex)
	if nVertex < 3:
	print("Vertex number is: ", nVertex)
	return
	else:
	ring = ogr.Geometry(ogr.wkbLinearRing)
	aCoords_gcs = np.full((nVertex,2), -9999.0, dtype=float)
	for j in range(nVertex):
	x1 = convert_360_to_180(aLonVertex[j])
	y1 = float(aLatVertex[j])
	ring.AddPoint(x1, y1)
	aCoords_gcs[j,0] = x1
	aCoords_gcs[j,1] = y1
	pass

	#add the closing point, which is the first vertex
	x1 = convert_360_to_180(aLonVertex[0])
	y1 = aLatVertex[0]
	ring.AddPoint(x1, y1) #double check
	pPolygon = ogr.Geometry(ogr.wkbPolygon)
	pPolygon.AddGeometry(ring)
	# Validate the geometry
	pPolygon_wkt = pPolygon.ExportToWkt()
	if not pPolygon.IsValid():
	print("Polygon is invalid...")
	#pPolygon = pPolygon.Buffer(0) # Buffering by 0 can fix some geometry issues
	#pPolygon = pPolygon.Simplify(0.0001) # Simplify with a tolerance
	else:
	print("Polygon is valid...")

	print("Polygon WKT:", pPolygon_wkt)
	if iFlag_plot_in is not None:
	#use a simple method to plot the polygon
	#get the path where this script is located
	sPath_script = os.path.dirname(os.path.abspath(__file__))
	sFilename_png = 'mesh_cell_' + str(lCellID_in) + '_culled.png'
	sFilename_png_out = os.path.join(sPath_script, sFilename_png)
	plot_wkt_polygon(pPolygon_wkt, sFilename_png_out)
	if sFilename_mpas_mesh_netcdf_base is not None:
	#find the base mesh cell using the center lon and lat
	lCell_index_base = np.where((aLongitudeCell_base == dLon_center_360) & (aLatitudeCell_base == dLat_center))
	#check we can find a base cell or not
	if len(lCell_index_base) == 1 and lCell_index_base[0] != 0:
	lCellID_base = aIndexToCellID_base[lCell_index_base]
	aVertexOnCellIndex = np.array(aVertexOnCell_base[lCell_index_base,:])
	dummy0 = np.where(aVertexOnCellIndex > 0)
	aVertexIndex = aVertexOnCellIndex[dummy0] - 1
	aLonVertex_base = aLongitudeVertex_base[aVertexIndex]
	aLatVertex_base = aLatitudeVertex_base[aVertexIndex]
	nVertex = len(aLonVertex)
	if nVertex < 3:
	print("Base vertex number is: ", nVertex)
	return
	else:
	ring = ogr.Geometry(ogr.wkbLinearRing)
	aCoords_gcs = np.full((nVertex,2), -9999.0, dtype=float)
	for j in range(nVertex):
	x1 = convert_360_to_180(aLonVertex_base[j])
	y1 = float(aLatVertex_base[j])
	ring.AddPoint(x1, y1)
	aCoords_gcs[j,0] = x1
	aCoords_gcs[j,1] = y1
	pass

	#add the closing point, which is the first vertex
	x1 = convert_360_to_180(aLonVertex_base[0])
	y1 = aLatVertex_base[0]
	ring.AddPoint(x1, y1) #double check
	pPolygon = ogr.Geometry(ogr.wkbPolygon)
	pPolygon.AddGeometry(ring)
	# Validate the geometry
	pPolygon_wkt_base = pPolygon.ExportToWkt()
	if not pPolygon.IsValid():
	print("Base polygon is invalid...")
	#pPolygon = pPolygon.Buffer(0) # Buffering by 0 can fix some geometry issues
	#pPolygon = pPolygon.Simplify(0.0001) # Simplify with a tolerance
	else:
	print("Base polygon is valid...")

	sFilename_png = 'mesh_cell_' + str(lCellID_in) + '_base.png'
	sFilename_png_out = os.path.join(sPath_script, sFilename_png)
	plot_wkt_polygon(pPolygon_wkt_base, sFilename_png_out)
	else:
	print("The base cell is not found...")
	continue
	pass



	return

	#create a main call
	if __name__ == '__main__':
	sFilename_mpas_mesh_netcdf_culled = '/compyfs/liao313/04model/pyflowline/conus/pyflowline20241201019/jigsaw/out/invert_mesh.nc'
	sFilename_mpas_mesh_netcdf_base = '/compyfs/liao313/04model/pyflowline/conus/pyflowline20241201019/jigsaw/out/base_mesh.nc'
	aCellID_in = list()
	#Polygon is invalid... 1672
	#Polygon WKT: POLYGON ((122.280664576308 18.2221473603637 0,122.262233339828 18.2160835887451 0,122.263243646392 18.1964473296003 0,122.279796144809 18.1912927994405 0,122.287660730398 18.2094542663177 0,122.278850506809 18.2232593941033 0,122.280664576308 18.2221473603637 0))
	#Polygon is invalid... 13994
	#Polygon WKT: POLYGON ((-63.9008296243521 50.3336863925212 0,-63.9273355026816 50.3213375919327 0,-63.9062274615777 50.3021870895725 0,-63.9271251813239 50.3195213928764 0,-63.8814145262508 50.3029476266043 0,-63.8690894249787 50.311717913351 0,-63.8704197044186 50.3233279141183 0,-63.9008296243521 50.3336863925212 0))
	#Polygon is invalid... 15428
	#Polygon WKT: POLYGON ((-48.7470719201964 61.4094980978841 0,-48.7078667093716 61.4115580717041 0,-48.7076968290549 61.4159935394359 0,-48.7244691494012 61.4238440398355 0,-48.7482711956022 61.4177357340897 0,-48.7500637709318 61.4113714878841 0,-48.7318223408571 61.407312966789 0,-48.7470719201964 61.4094980978841 0))
	#Polygon is invalid... 17062
	#Polygon WKT: POLYGON ((-129.539599690119 70.1323831312393 0,-129.535795196565 70.124261028838 0,-129.512483903554 70.1172844555891 0,-129.465730846556 70.1172188707654 0,-129.459111064388 70.1364237439907 0,-129.465652367021 70.1172520619199 0,-129.498714599157 70.149964048113 0,-129.539599690119 70.1323831312393 0))
	#Polygon is invalid... 17793
	#Polygon WKT: POLYGON ((-71.4446320396707 10.4020008778489 0,-71.4516631039428 10.3961118391854 0,-71.456830793219 10.4214327874876 0,-71.4760303796148 10.4223603510441 0,-71.4842067505045 10.4056232943428 0,-71.4713271460859 10.3935027872622 0,-71.4446320396707 10.4020008778489 0))
	#Polygon is invalid... 24092
	#Polygon WKT: POLYGON ((-15.5382028906378 11.5876827405335 0,-15.535595747457 11.5850973512438 0,-15.5099707980737 11.591618233233 0,-15.5071308491261 11.599961078997 0,-15.5246821263325 11.613672338718 0,-15.5098012098433 11.6054745495611 0,-15.5420411080539 11.6092003139624 0,-15.5382028906378 11.5876827405335 0))
	#Polygon is invalid... 24513
	#Polygon WKT: POLYGON ((-43.3410601521804 -2.35730294671818 0,-43.3101695695966 -2.36488149251201 0,-43.3081346848379 -2.35641950002847 0,-43.3166767565955 -2.33487861443757 0,-43.3342404233832 -2.3423723030464 0,-43.3405705233804 -2.36152207379487 0,-43.3410601521804 -2.35730294671818 0))
	#Polygon is invalid... 33144
	#Polygon WKT: POLYGON ((119.364514438706 -1.17964356953732 0,119.348411574817 -1.17792521472226 0,119.339621749198 -1.18809595090908 0,119.34848087937 -1.20993515223943 0,119.34738302738 -1.2089944663736 0,119.365480682366 -1.20127151411875 0,119.364514438706 -1.17964356953732 0))
	#Polygon is invalid... 35167
	#Polygon WKT: POLYGON ((-35.9944101244893 -5.05398122354815 0,-35.9801257279569 -5.06874982595671 0,-35.9656043467338 -5.06703324195992 0,-35.9708716349697 -5.04776250575661 0,-35.9648702774745 -5.06163712690284 0,-35.9884656481591 -5.03975398282233 0,-35.9944101244893 -5.05398122354815 0))
	#Polygon is invalid... 35616
	#Polygon WKT: POLYGON ((-157.006243393183 71.2109730489155 0,-157.010511724013 71.1907258531127 0,-156.984279984155 71.1778101013443 0,-156.986514943451 71.1781925881701 0,-156.93530954857 71.1846969871186 0,-156.948468763097 71.2107407300836 0,-157.006243393183 71.2109730489155 0))
	#Polygon is invalid... 35641
	#Polygon WKT: POLYGON ((33.7201130735598 28.0767992360369 0,33.7389806225413 28.0809644371193 0,33.734809560818 28.0775273580535 0,33.734179155333 28.0972618198086 0,33.7132054411147 28.1016558473581 0,33.7023215895678 28.086756860363 0,33.7201130735598 28.0767992360369 0))
	#Polygon is invalid... 38117
	#Polygon WKT: POLYGON ((148.343510865271 -5.59033620753651 0,148.355315036408 -5.57459782311827 0,148.3452916253 -5.55419706603241 0,148.343320630905 -5.55389792166596 0,148.329597677703 -5.5616335416741 0,148.326235132406 -5.57660836257517 0,148.345384562197 -5.59062056647374 0,148.343510865271 -5.59033620753651 0))
	#Polygon is invalid... 41549
	#Polygon WKT: POLYGON ((137.74958112273 -35.1129190969386 0,137.728498752014 -35.1157816659266 0,137.763899893261 -35.0954158510455 0,137.735133286168 -35.0906702014205 0,137.719432127827 -35.1031570555177 0,137.722144962911 -35.1135022577956 0,137.74958112273 -35.1129190969386 0))
	#Polygon is invalid... 41660
	#Polygon WKT: POLYGON ((-16.1689593883591 12.2942482599369 0,-16.1656619757907 12.2913379381728 0,-16.1501863914894 12.2909515078012 0,-16.1388797897555 12.302421020704 0,-16.1457460487068 12.3239541811496 0,-16.1435360602088 12.3216513216158 0,-16.1733726860637 12.3104519598962 0,-16.1689593883591 12.2942482599369 0))
	#Polygon is invalid... 42129
	#Polygon WKT: POLYGON ((-71.873745052447 10.2041331098005 0,-71.8695608031783 10.2010706472204 0,-71.8505573208394 10.2005214388756 0,-71.8431490329504 10.2148235313675 0,-71.8501723767827 10.2314272143016 0,-71.8495997694914 10.2310373849819 0,-71.8720248298625 10.2204736226489 0,-71.873745052447 10.2041331098005 0))
	#Polygon is invalid... 45022
	#Polygon WKT: POLYGON ((-3.49990419525511 50.4458670833471 0,-3.50109363562478 50.4650696303941 0,-3.53783879506921 50.4773819921189 0,-3.52408079611428 50.4765892264093 0,-3.54618447523524 50.4521880445044 0,-3.52947657697973 50.4430455316903 0,-3.49990419525511 50.4458670833471 0))
	#Polygon is invalid... 46969
	#Polygon WKT: POLYGON ((-42.8889909455578 61.0765107756435 0,-42.9121083332119 61.0772667813284 0,-42.916511986311 61.0595397260809 0,-42.8723073974284 61.0505799864484 0,-42.869995504361 61.0608183805887 0,-42.8707395979006 61.0514417996377 0,-42.8889909455578 61.0765107756435 0))
	#Polygon is invalid... 50195
	#Polygon WKT: POLYGON ((-94.368666425969 68.9908495945327 0,-94.3175518357232 68.9953417755255 0,-94.3815668165992 68.9736632134828 0,-94.3425952667274 68.9661781569792 0,-94.2997990225276 68.9815048216357 0,-94.313890214263 68.9949962318141 0,-94.368666425969 68.9908495945327 0))
	#Polygon is invalid... 51427
	#Polygon WKT: POLYGON ((-91.9305657042177 29.8364298178733 0,-91.9471611860289 29.8502110560581 0,-91.9683012276131 29.8435993051381 0,-91.9705116516197 29.826359785042 0,-91.966002264764 29.8170637700613 0,-91.9592633280558 29.8132609164735 0,-91.931947199102 29.8332217173152 0,-91.9315475545814 29.8412767907539 0,-91.9305657042177 29.8364298178733 0))
	#Polygon is invalid... 52022
	#Polygon WKT: POLYGON ((119.452814084782 26.0049533277492 0,119.461850081725 25.9959711700541 0,119.484939468505 25.9981496514189 0,119.490311905252 26.0107563312908 0,119.475593022448 26.0293297088006 0,119.481284836318 26.0265138358557 0,119.452814084782 26.0049533277492 0))
	#Polygon is invalid... 53684
	#Polygon WKT: POLYGON ((134.880996223007 34.6996946747361 0,134.885760101184 34.7052291283607 0,134.86918556633 34.7273923580406 0,134.846823160703 34.7248603694141 0,134.846104671205 34.7037126681144 0,134.84636206586 34.724189250281 0,134.861751304972 34.6951355631687 0,134.880996223007 34.6996946747361 0))
	#Polygon is invalid... 53873
	#Polygon WKT: POLYGON ((67.8621346406756 71.4536132089809 0,67.8603493762603 71.4692951262031 0,67.820678200099 71.4835092147059 0,67.7810267313965 71.4725449644734 0,67.7937258331278 71.4477814842135 0,67.7825844661424 71.4529756746679 0,67.8476571359303 71.4481797847844 0,67.8621346406756 71.4536132089809 0))
	#Polygon is invalid... 54221
	#Polygon WKT: POLYGON ((-42.3203661613626 62.709446071854 0,-42.3313270148944 62.7199910936606 0,-42.3741515666778 62.7235736362693 0,-42.3760364767373 62.7051599709813 0,-42.3803038139625 62.7147640854633 0,-42.3462106371608 62.6948465209412 0,-42.3274484550065 62.6974134802222 0,-42.3203661613626 62.709446071854 0))
	#Polygon is invalid... 55261
	#Polygon WKT: POLYGON ((-96.2548737107388 68.5209704407071 0,-96.2977960367591 68.5145042550792 0,-96.3095326985719 68.4931769995066 0,-96.2957756974138 68.4880380626558 0,-96.2514869022441 68.4879607107365 0,-96.2306622412806 68.5028910535988 0,-96.2581634506687 68.5217217190116 0,-96.2548737107388 68.5209704407071 0))
	#Polygon is invalid... 57712
	#aCellID_in.append(1672+1)
	#aCellID_in.append(13994+1)
	#aCellID_in.append(15428+1)
	aCellID_in.append(17062+1)
	#aCellID_in.append(17793+1)
	#aCellID_in.append(24092+1)
	#aCellID_in.append(24513+1)
	#aCellID_in.append(33144+1)
	#aCellID_in.append(35167+1)

	check_mesh_cell_polygon(sFilename_mpas_mesh_netcdf_culled, aCellID_in, iFlag_plot_in = 1, sFilename_mpas_mesh_netcdf_base = sFilename_mpas_mesh_netcdf_base )