Try to write FVCOM mesh as geojson



In [1]:

    
from shapely.geometry import MultiPolygon, mapping, polygon
import json



In [2]:

    
#url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc'
#url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_FVCOM_OCEAN_MASSBAY_FORECAST.nc'
url = 'http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_WAVE_FORECAST.nc'

use netcdf4 because UGRID takes longer

import pyugrid ug = pyugrid.UGrid.from_ncfile(url) lon = ug.nodes[:,0] lat = ug.nodes[:,1] nv = ug.faces[:]



In [3]:

    
import netCDF4
nc = netCDF4.Dataset(url)
ncv = nc.variables
lon = ncv['lon'][:]
lat = ncv['lat'][:]
nv = ncv['nv'][:,:].T - 1



In [4]:

    
#mp = MultiPolygon([polygon.Polygon(zip(lon[element],lat[element])) for element in nv])
mp = MultiPolygon([polygon.Polygon(zip(lon[element],lat[element])) for element in nv[0:5]])



In [5]:

    
type(mp)









    Out[5]:





shapely.geometry.multipolygon.MultiPolygon



In [6]:

    
mp









    



WARNING:shapely.geos:Self-intersection at or near point -59.825973510742188 46.058303833007812






    Out[6]:



In [7]:

    
json.dumps(mapping(mp))









    Out[7]:





'{"type": "MultiPolygon", "coordinates": [[[[-59.890201568603516, 46.09441375732422], [-59.81068801879883, 46.14595031738281], [-59.82597351074219, 46.05830383300781], [-59.890201568603516, 46.09441375732422]]], [[[-59.82597351074219, 46.05830383300781], [-59.81068801879883, 46.14595031738281], [-59.75267791748047, 46.090110778808594], [-59.82597351074219, 46.05830383300781]]], [[[-59.701324462890625, 46.038551330566406], [-59.82597351074219, 46.05830383300781], [-59.75267791748047, 46.090110778808594], [-59.701324462890625, 46.038551330566406]]], [[[-59.82597351074219, 46.05830383300781], [-59.701324462890625, 46.038551330566406], [-59.788787841796875, 46.0040283203125], [-59.82597351074219, 46.05830383300781]]], [[[-59.788787841796875, 46.0040283203125], [-59.701324462890625, 46.038551330566406], [-59.729705810546875, 45.940338134765625], [-59.788787841796875, 46.0040283203125]]]]}'



In [8]:

    
with open('ugrid.json','w') as f:
    json.dump(mapping(mp), f)



In [8]:



In [8]: