Math Forum Discussions

rachael

Posts: 16
Registered: 2/10/09

Re: griddata error: qhull precision error, ERRONEOUS FACET
Posted: Feb 24, 2009 3:52 PM

Plain Text

"rachael " <rmueller@coas.oregonstate.edu> wrote in message <gnso7i$k98$1@fred.mathworks.com>...
> I am trying to grid data 1 arc second grid onto a 4 Km grid, using griddata. I have tried inputting the data twice, and I get a similar error both times. The first time, I used a version of the data the comes from someone else; the second time, I am using a version of the data that I created from a downloaded dataset. The data that I "created" is from the GEBCO dataset of gridded bathymetry.
> <http://www.gebco.net/data_and_products/gridded_bathymetry_data/>
> It's 1 arc-minute, with equi-distant grid spacing. I create my latitude (LAT), longitude(LON) and bathymetry (Z) grids as follows:
>
> > ncload('/Users/rmueller/Data/GEBCO /62s75s80w40w_GEBCO_1min.nc','x_range','y_range','z_range','z');
>
> > ncload('/Users/rmueller/Data/GEBCO/62s75s80w40w_GEBCO_1min.nc','spacing','dimension');
>
> > dx = spacing(1);dy=spacing(2);
>
> > x=[x_range(2):-1*dx:x_range(1)];y=[y_range(2):-1*dy:y_range(1)];
>
> > [LAT,LON]=meshgrid(y,x);
>
> > Z = flipud(reshape(z,dimension(1),dimension(2)));
>
> The errors arise when I use griddata to map this data onto a 4km spaced, polar stereographic, lat_rho/lon_rho grid using the following
> >z = griddata(LON,LAT,Z,lon_rho,lat_rho);
>
> I have used this lat_rho, lon_rho to convert a different dataset w/o problems, so I think it's the GEBCO data. Since I get a similar error with the data as converted into matlab format by two different users, I think the problem is more due to the data than the format I have it in matlab. But what is the problem??? Is the 1 arc-minute too great of a grid spacing to map onto a 4 Km grid???
>
> The two errors are shown in detail below.
>
> ERROR 1:
> ??? Error using ==> qhullmx
> qhull precision error: initial facet 2 is coplanar with the interior point
>
> ERRONEOUS FACET:
> While executing: | qhull d Qt Qbb Qc
> Options selected for Qhull 2003.1 2003/12/30:
> delaunay Qtriangulate Qbbound-last Qcoplanar-keep _pre-merge
> _zero-centrum Pgood Qinterior-keep _max-width 35 Error-roundoff 1.1e-13
> _one-merge 7.4e-13 Visible-distance 2.1e-13 U-coplanar-distance 2.1e-13
> Width-outside 4.2e-13 _wide-facet 1.3e-12
> precision problems (corrected unless 'Q0' or an error)
> 2
> flipped facets
> The input to qhull appears to be less than 3 dimensional, or a
> computation has overflowed.
> Qhull could not construct a clearly convex simplex from points:
> The center point is coplanar with a facet, or a vertex is coplanar
> with a neighboring facet. The maximum round off error for
> computing distances is 1.1e-13. The center point, facets and distances
> to the center point are as follows:
>
> facet
> p2016960
> p2100
> p0
> distance= -3.6e-15
>
> facet
> p2019060
> p2100
> p0
> distance= 0
>
> facet
> p2019060
> p2016960
> p0
> distance= -1.4e-14
>
> facet
> p2019060
> p2016960
> p2100
> distance= -5.3e-15
>
> These points either have a maximum or minimum x-coordinate, or
> they maximize the determinant for k coordinates. Trial points
> are first selected from points that maximize a coordinate.
> The min and max coordinates for each dimension are:
> 0: -75 -2.225e-308 difference= 75
> 1: -76 -2.225e-308 difference= 76
> 2: -3.553e-15 35 difference= 35
> If the input should be full dimensional, you have several options that
> may determine an initial simplex:
> - use 'QJ' to joggle the input and make it full dimensional
> - use 'QbB' to scale the points to the unit cube
> - use 'QR0' to randomly rotate the input for different maximum points
> - use 'Qs' to search all points for the initial simplex
> - use 'En' to specify a maximum roundoff error less than 1.1e-13.
> - trace execution with 'T3' to see the determinant for each point.
> If the input is lower dimensional:
> - use 'QJ' to joggle the input and make it full dimensional
> - use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
> pick the coordinate with the least range. The hull will have the
> correct topology.
> - determine the flat containing the points, rotate the points
> into a coordinate plane, and delete the other coordinates.
> - add one or more points to make the input full dimensional.
> This is a Delaunay triangulation and the input is co-circular or co-spherical:
> - use 'Qz' to add a point "at infinity" (i.e., above the paraboloid)
> - or use 'QJ' to joggle the input and avoid co-circular data
> Error in ==> delaunayn at 117
> t = qhullmx(x', 'd ', opt);
> Error in ==> griddata>linear at 151
> tri = delaunayn([x y]);
> Error in ==> griddata at 120
> zi = linear(x,y,z,xi,yi,opt);
>
>
> ERROR 2:
> ??? Error using ==> qhullmx
> While executing: | qhull d Qt Qbb Qc
> Options selected for Qhull 2003.1 2003/12/30:
> delaunay Qtriangulate Qbbound-last Qcoplanar-keep _pre-merge
> _zero-centrum Pgood Qinterior-keep _max-width 40 Error-roundoff 1.1e-13
> _one-merge 7.8e-13 Visible-distance 2.2e-13 U-coplanar-distance 2.2e-13
> Width-outside 4.4e-13 _wide-facet 1.3e-12
> Last point added to hull was p383250.
> Last merge was #3109857.
> Qhull has finished constructing the hull.
> At error exit:
> precision problems (corrected unless 'Q0' or an error)
> 3109857
> coplanar horizon facets for new vertices
> Error in ==> delaunayn at 117
> t = qhullmx(x', 'd ', opt);
> Error in ==> griddata>linear at 151
> tri = delaunayn([x y]);
> Error in ==> griddata at 120
> zi = linear(x,y,z,xi,yi,opt);
>
> Thanks.

Problem solved! Now that I know the answer, the cryptic error message makes more sense. The 'QJ' mentioned above is used as an "option" in griddata. The solution that worked for me is:

z = griddata(lon_g,lat_g,z_g,lon_rho,lat_rho, 'linear',{'QJ'});