¿Alisar polígonos en el mapa de contorno?

52

Aquí hay un mapa de contorno para el que están disponibles todos los polígonos de niveles.

Vamos a preguntar cómo suavizar los polígonos manteniendo todos los vértices preservados en sus ubicaciones exactas.

De hecho, el contorno se realiza sobre los datos de una cuadrícula, puede sugerir que luego suavice los datos de la cuadrícula y, por lo tanto, el contorno resultante será más suave. Tenga en cuenta que esto no funciona como lo deseo, ya que la función de suavizado como el filtro Gaussiano eliminará pequeños paquetes de datos y cambiará el rango de la tercera variable, por ejemplo, la altura que no está permitida en mi aplicación.

En realidad, estoy buscando un fragmento de código (preferiblemente en Python ) que pueda suavizar los polígonos 2D (cualquier tipo: convexo, cóncavo, auto-intersectorial, etc.) razonablemente indoloro (olvide las páginas de códigos) y preciso.

Para su información, hay una función en ArcGIS que hace esto perfectamente, pero el uso de aplicaciones comerciales de terceros no es mi elección para esta pregunta.

ingrese la descripción de la imagen aquí


1)

Scipy.interpolate:

ingrese la descripción de la imagen aquí

Como puede ver, las splines resultantes (rojo) no son satisfactorias.

2)

Aquí está el resultado usando el código dado aquí . ¡No está funcionando bien!

ingrese la descripción de la imagen aquí

3)

Para mí, la mejor solución debería ser algo como la siguiente figura en la que un cuadrado se suaviza gradualmente cambiando solo un valor. Espero un concepto similar para suavizar cualquier forma de polígonos.

ingrese la descripción de la imagen aquí

Satisfaciendo la condición de que la spline pase los puntos:

ingrese la descripción de la imagen aquí

4)

Aquí está mi implementación de "idea de whuber" línea por línea en Python en sus datos. Posiblemente hay algunos errores ya que los resultados no son buenos.

ingrese la descripción de la imagen aquí

K = 2 es un desastre y por lo tanto para k> = 4.

5)

Eliminé un punto en la ubicación problemática y la spline resultante ahora es idéntica a la de Whuber. Pero todavía es una pregunta de por qué el método no funciona en todos los casos.

ingrese la descripción de la imagen aquí

6)

Un buen suavizado para los datos de Whuber puede ser el siguiente (dibujado por un software de gráficos vectoriales) en el que se ha agregado un punto extra sin problemas (compárese con la actualización

4):

ingrese la descripción de la imagen aquí

7)

Vea el resultado de la versión de Python del código de Whuber para algunas formas icónicas:

ingrese la descripción de la imagen aquí
Tenga en cuenta que el método parece no funcionar para las polilíneas. Para la esquina, la polilínea (contorno) verde es lo que quiero pero se puso rojo. Esto debe abordarse ya que los mapas de contorno son siempre polilíneas, aunque las polilíneas cerradas pueden tratarse como polígonos como en mis ejemplos. Tampoco es que el problema surgido en la actualización 4 aún no se haya abordado.

8) [mi último]

Aquí está la solución final (¡no perfecta!):

ingrese la descripción de la imagen aquí

Recuerde que tendrá que hacer algo sobre el área señalada por las estrellas. Quizás haya un error en mi código o el método propuesto necesita un mayor desarrollo para considerar todas las situaciones y proporcionar los resultados deseados.

Desarrollador
fuente
¿Cómo estás generando contornos 'polígonos'? ¿No serían siempre líneas, ya que un contorno que se cruza con el borde de un DEM nunca se cerraría sobre sí mismo?
pistacho
He utilizado la función v.generalize en GRASS para suavizar las líneas de contorno con resultados decentes, aunque puede tomar un tiempo para mapas con contornos muy densos.
pistacho
@pistachionut Puede considerar que los niveles de contorno son polilíneas. Estoy buscando código puro en la primera etapa. Si no está disponible, paquete ligero para Python.
Desarrollador
Quizás mire scipy.org/Cookbook/Interpolation porque suena como si quisiera spline
PolyGeo
1
La curva @Pablo Bezier en su enlace funciona bien para polilíneas. Whuber's funciona casi bien para polígonos. Para que juntos puedan abordar la pregunta. Muchas gracias por compartir tus conocimientos de forma gratuita.
Desarrollador

Respuestas:

37

La mayoría de los métodos para spline secuencias de números spline polígonos. El truco es hacer que las splines se "cierren" suavemente en los puntos finales. Para hacer esto, "envuelve" los vértices alrededor de los extremos. Luego spline las coordenadas x e y por separado.

Aquí hay un ejemplo de trabajo en R. Utiliza el splineprocedimiento cúbico predeterminado disponible en el paquete de estadísticas básicas. Para obtener más control, sustituya casi cualquier procedimiento que prefiera: solo asegúrese de que se divide entre los números (es decir, los interpola) en lugar de simplemente usarlos como "puntos de control".

#
# Splining a polygon.
#
#   The rows of 'xy' give coordinates of the boundary vertices, in order.
#   'vertices' is the number of spline vertices to create.
#              (Not all are used: some are clipped from the ends.)
#   'k' is the number of points to wrap around the ends to obtain
#       a smooth periodic spline.
#
#   Returns an array of points. 
# 
spline.poly <- function(xy, vertices, k=3, ...) {
    # Assert: xy is an n by 2 matrix with n >= k.

    # Wrap k vertices around each end.
    n <- dim(xy)[1]
    if (k >= 1) {
        data <- rbind(xy[(n-k+1):n,], xy, xy[1:k, ])
    } else {
        data <- xy
    }

    # Spline the x and y coordinates.
    data.spline <- spline(1:(n+2*k), data[,1], n=vertices, ...)
    x <- data.spline$x
    x1 <- data.spline$y
    x2 <- spline(1:(n+2*k), data[,2], n=vertices, ...)$y

    # Retain only the middle part.
    cbind(x1, x2)[k < x & x <= n+k, ]
}

Para ilustrar su uso, creemos un polígono pequeño (pero complicado).

#
# Example polygon, randomly generated.
#
set.seed(17)
n.vertices <- 10
theta <- (runif(n.vertices) + 1:n.vertices - 1) * 2 * pi / n.vertices
r <- rgamma(n.vertices, shape=3)
xy <- cbind(cos(theta) * r, sin(theta) * r)

Spline usando el código anterior. Para suavizar la spline, aumente el número de vértices de 100; para que sea menos suave, disminuya el número de vértices.

s <- spline.poly(xy, 100, k=3)

Para ver los resultados, trazamos (a) el polígono original en rojo discontinuo, mostrando el espacio entre el primer y el último vértice (es decir, sin cerrar su polilínea límite); y (b) la spline en gris, mostrando una vez más su espacio. (Debido a que el espacio es tan pequeño, sus puntos finales se resaltan con puntos azules).

plot(s, type="l", lwd=2, col="Gray")
lines(xy, col="Red", lty=2, lwd=2)
points(xy, col="Red", pch=19)
points(s, cex=0.8)
points(s[c(1,dim(s)[1]),], col="Blue", pch=19)

Polígono estriado

whuber
fuente
55
Buena respuesta. ¿Hay alguna forma de garantizar que los contornos no se crucen como resultado del suavizado?
Kirk Kuykendall
Esa es una buena pregunta, @Kirk. No conozco ningún método para garantizar el no cruce de esta forma de suavizado. (De hecho, ni siquiera veo cómo garantizar que la polilínea suavizada no se cruce entre sí. Sin embargo, esto no es un gran problema para la mayoría de los contornos). Para hacer eso, necesitaría volver al original DEM y en su lugar use un mejor método para calcular los contornos en primer lugar. (No son mejores métodos - que se conocen desde hace mucho tiempo -. Que yo sepa, pero algunos de los más populares GISes no usarlos)
whuber
Primero, todavía estoy trabajando para implementar su respuesta en Python, pero no he tenido éxito. Segundo, ¿cuál será el resultado si aplica su método en un cuadrado? Puede referirse a los que he dibujado en la pregunta.
Desarrollador
1
Acepté esto como respuesta ya que da una buena solución. Aunque no es perfecto pero me dio algunas ideas para trabajar, espero encontrar una solución que satisfaga los puntos que mencioné anteriormente en mi pregunta y comentarios. También puede considerar los comentarios de Whuber para la pregunta [QC], hay buenos trucos allí. Por último, debo decir que la traducción a Python es casi sencilla, ya que tiene instalado el encantador paquete Scipy. Considere también el comentario de Pablo en QC como posible solución para polilíneas, es decir, curvas de Bezier. Buena suerte a todos.
Desarrollador
1
Al ver sus respuestas, lamento no haberme ocupado bien de mis matemáticas.
vinayan
2

Sé que esta es una publicación antigua, pero apareció en Google por algo que estaba buscando, así que pensé en publicar mi solución.

No veo esto como un ejercicio de ajuste de curva 2D, sino más bien como uno 3D. Al considerar los datos como 3D, podemos asegurarnos de que las curvas nunca se crucen entre sí, y podemos usar la información de otros contornos para mejorar nuestra estimación para la actual.

El siguiente extracto de iPython utiliza la interpolación cúbica proporcionada por SciPy. Tenga en cuenta que los valores z que he trazado no son importantes, siempre y cuando todos los contornos sean equidistantes en altura.

In [1]: %pylab inline
        pylab.rcParams['figure.figsize'] = (10, 10)
        Populating the interactive namespace from numpy and matplotlib

In [2]: import scipy.interpolate as si

        xs = np.array([0.0, 0.0, 4.5, 4.5,
                       0.3, 1.5, 2.3, 3.8, 3.7, 2.3,
                       1.5, 2.2, 2.8, 2.2,
                       2.1, 2.2, 2.3])
        ys = np.array([0.0, 3.0, 3.0, 0.0,
                       1.1, 2.3, 2.5, 2.3, 1.1, 0.5,
                       1.1, 2.1, 1.1, 0.8,
                       1.1, 1.3, 1.1])
        zs = np.array([0,   0,   0,   0,
                       1,   1,   1,   1,   1,   1,
                       2,   2,   2,   2,
                       3,   3,   3])
        pts = np.array([xs, ys]).transpose()

        # set up a grid for us to resample onto
        nx, ny = (100, 100)
        xrange = np.linspace(np.min(xs[zs!=0])-0.1, np.max(xs[zs!=0])+0.1, nx)
        yrange = np.linspace(np.min(ys[zs!=0])-0.1, np.max(ys[zs!=0])+0.1, ny)
        xv, yv = np.meshgrid(xrange, yrange)
        ptv = np.array([xv, yv]).transpose()

        # interpolate over the grid
        out = si.griddata(pts, zs, ptv, method='cubic').transpose()

        def close(vals):
            return np.concatenate((vals, [vals[0]]))

        # plot the results
        levels = [1, 2, 3]
        plt.plot(close(xs[zs==1]), close(ys[zs==1]))
        plt.plot(close(xs[zs==2]), close(ys[zs==2]))
        plt.plot(close(xs[zs==3]), close(ys[zs==3]))
        plt.contour(xrange, yrange, out, levels)
        plt.show()

Resultado cúbico interpolado

Los resultados aquí no se ven los mejores, pero con tan pocos puntos de control siguen siendo perfectamente válidos. Observe cómo se estira la línea verde ajustada para seguir el contorno azul más ancho.

Gilly
fuente
Las curvas suaves ajustadas deben permanecer lo más cerca posible del polígono / polilínea original.
Desarrollador
1

Escribí casi exactamente el paquete que está buscando ... pero estaba en Perl, y fue hace más de una década: GD :: Polyline . Utilizaba curvas de Bezier cúbicas 2D y "suavizaba" un polígono o "polilínea" arbitrario (mi nombre para lo que ahora se llama comúnmente "LineString").

El algoritmo consistió en dos pasos: dados los puntos en el Polígono, agregue dos puntos de control Bezier entre cada punto; luego llame a un algoritmo simple para hacer una aproximación por partes de la spline.

La segunda parte es fácil; La primera parte fue un poco de arte. Aquí fue la idea: considerar un "segmento de control de" un Vertex N: vN. El segmento de control fue de tres puntos colineales: [cNa, vN, cNb]. El punto central era el vértice. La pendiente de este control seg fue igual a la pendiente del Vértice N-1 al Vértice N + 1. La longitud de la porción izquierda de este segmento era 1/3 de la longitud del Vértice N-1 al Vértex N, y la longitud de la porción derecha de este segmento era 1/3 de la longitud del Vértice N al Vértex N + 1.

Si la curva original era cuatro vértices: [v1, v2, v3, v4]entonces el cada vértice ahora conseguir un segmento de control de la forma: [c2a, v2, c2b]. Une estos juntos de esta manera: [v1, c1b, c2a, v2, c2b, c3a, v3, c3b, c4a, v4]y mójalos de cuatro en cuatro como los cuatro puntos de Bezier:, [v1, c1b, c2a, v2]luego [v2, c2b, c3a, v3], y así sucesivamente. Como [c2a, v2, c2b]eran co-lineales, la curva resultante será suave en cada vértice.

Por lo tanto, esto también cumple con su requisito de parametrizar la "tensión" de la curva: use un valor menor que 1/3 para una curva "más ajustada", una más grande para un ajuste "más circular". En cualquier caso, la curva resultante siempre pasa por los puntos dados originales.

Esto dio como resultado una curva suave que "circunscribió" el Polígono original. También tenía alguna forma de "inscribir" una curva suave ... pero no veo eso en el código CPAN.

De todos modos, en este momento no tengo una versión disponible en Python, ni tengo ninguna figura. PERO ... si / cuando porto esto a Python, me aseguraré de publicar aquí.

Dan H
fuente
No se puede evaluar el código Perl, si es posible, agregue gráficos para demostrar cómo funcionaba.
Desarrollador