Ver Tema

lis · 2008-05-27 21:54 @954

Hola, mi problema es el siguiente: tengo una matriz (con extensión .dat) de 49x50, que en este caso son datos de temperatura, y debo hacer una interpolación bilineal a cualquier punto (x,y), en donde me muestre los 4 puntos que rodean a (x,y) y además me muestre el resultado de la interpolación.

Si alguien me puede ayudar con este script en Perl se lo agradecería mucho porque llevo en esto casi 3 semanas.

Muchas gracias...

explorer · 2008-05-28 06:16 @303

¿Quieres trabajar con toda la imagen de golpe o solo un pixel que elijas tu? ¿Cómo quieres elegirlo, con el ratón o con el teclado? ¿Y la salida, cómo debe ser?

Actualización: Ya encontré la expresión matemática de la interpolación bilineal.

explorer · 2008-05-28 08:03 @377

Este es un ejemplo de cálculo de interpolación bilineal:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

#!/usr/bin/perl

use warnings;

use strict;

# Toma de datos.

# Suponemos que el origen de coordenadas está en

# la coordenada superior izquierda de la tabla.

my $temperaturas = [

    [ 0, 0, 0, 0, 0, 0, 0, ],

    [ 0, 5, 4, 4, 6, 5, 0, ],

    [ 1, 4, 5, 5, 7, 6, 0, ],

    [ 3, 3, 4, 5, 4, 4, 0, ],

    [ 1, 2, 2, 3, 2, 1, 0, ],

    [ 0, 1, 0, 1, 1, 0, 0, ],

    [ 0, 0, 0, 0, 0, 0, 0, ],

    [ 0, 0, 0, 0, 0, 0, 0, ],

];

# Dimensiones del campo muestreado

my $alto_tabla_temperaturas  = @{$temperaturas     }-1; # Quitamos 1 para evitar el caso

my $ancho_tabla_temperaturas = @{$temperaturas->[0]}-1; # del fallo en los bordes

# Dimensiones físicas del campo

my $factor_ancho = 10;

my $factor_alto  = 15;

my $ancho = $ancho_tabla_temperaturas * $factor_ancho;

my $alto  = $alto_tabla_temperaturas  * $factor_alto;

# Cálculo de la temperatura interpolada bilinealmente

sub temperatura_interpolada {

    my ( $x, $y ) = @_;         # Posición dentro del campo físico

    return 'ERROR: posición fuera de campo'

        if $x < 0 or $y < 0 or $x >= $ancho or $y >= $alto;

    # Cálculo de las coordenadas de los píxeles que rodean

    # al indicado, dentro de la tabla de muestreo

    my ($a , $b ) = ( $x / $factor_ancho, $y / $factor_alto );

    my ($x1, $y1) = ( int $a            , int $b            );

    my ($x2, $y2) = ( $x1 + 1           , $y1 + 1           );

    # Cálculo de la posición de entrada, referida dentro del cuadro

    my ($xp, $yp) = ( $a - $x1, $b - $y1 );

    # Extraemos los valores muestreados del cuadro que lo rodea

    my ( $f00, $f10, $f01, $f11 ) = (

        $temperaturas->[$y1][$x1],

        $temperaturas->[$y1][$x2],

        $temperaturas->[$y2][$x1],

        $temperaturas->[$y2][$x2],

    );

    print "Coordenadas y valores de los puntos que le rodean, en la tabla de muestreo:\n";

    print "\t[ $x1, $y1 ]:$f00 --- [ $x2, $y1 ]:$f10\n";

    print "\t[ $x1, $y2 ]:$f01 --- [ $x2, $y2 ]:$f11\n";

    print "Posición relativa dentro del cuadro: $xp,$yp\n";

    # Cálculo del valor interpolado

    my ( $b1, $b2, $b3, $b4 ) = ( $f00, $f10-$f00, $f01-$f00, $f00-$f10-$f01+$f11 );

    return $b1  +  $b2 * $xp  +  $b3 * $yp  +  $b4 * $xp * $yp;

}

# Bucle principal

print "Cálculo de temperaturas\n";

print "Dimensiones del campo: $ancho x $alto\n";

print "Dimensiones de la muestra: $ancho_tabla_temperaturas x $alto_tabla_temperaturas\n";

print "Factor: $factor_ancho, $factor_alto\n";

while ( 1 ) {

    print 'Introduzca las coordenadas a calcular (x,y): ';

    my $entrada = <>;

    chomp $entrada;

    my ($x,$y) = split q{,}, $entrada;

    $x += 0; $y += 0;

    next if !$x or !$y;

    my $t = temperatura_interpolada($x,$y);

    print "Cálculo del valor interpolado: $t\n";

    print "---------------\n";

}

__END__

Una posible salida sería:

Código:

Cálculo de temperaturas
Dimensiones del campo: 60 x 105
Dimensiones de la muestra: 6 x 7
Factor: 10, 15
Introduzca las coordenadas a calcular (x,y): 37,14
Coordenadas y valores de los puntos que le rodean, en la tabla de muestreo:
        [ 3, 0 ]:0 --- [ 4, 0 ]:0
        [ 3, 1 ]:4 --- [ 4, 1 ]:6
Posición relativa dentro del cuadro: 0.7,0.933333333333333
Cálculo del valor interpolado: 5.04

explorer · 2008-05-28 11:29 @520

Obviamente, es mucho mejor usar PDL, para hacer cosas serias...

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

#!/usr/bin/perl

use warnings;

use strict;

use PDL;

use PDL::Image2D;

my $temperaturas = pdl([

    [ 0, 0, 0, 0, 0, 0, 0, ],

    [ 0, 5, 4, 4, 6, 5, 0, ],

    [ 1, 4, 5, 5, 7, 6, 0, ],

    [ 3, 3, 4, 5, 4, 4, 0, ],

    [ 1, 2, 2, 3, 2, 1, 0, ],

    [ 0, 1, 0, 1, 1, 0, 0, ],

    [ 0, 0, 0, 0, 0, 0, 0, ],

    [ 0, 0, 0, 0, 0, 0, 0, ],

]);

print $temperaturas;

my $superficie_real = zeroes(14,14);

bilin2d($temperaturas, $superficie_real);

print $superficie_real;

__END__

Código:

[
 [0 0 0 0 0 0 0]
 [0 5 4 4 6 5 0]
 [1 4 5 5 7 6 0]
 [3 3 4 5 4 4 0]
 [1 2 2 3 2 1 0]
 [0 1 0 1 1 0 0]
 [0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0]
]

[
 [             0              0              0              0              0              0              0              0              0              0              0              0              0              0]
 [             0      1.2426036      2.4852071      2.4852071      2.2366864      2.1538462      2.1538462      2.4023669      2.8994083       3.147929      2.8994083      2.4852071      1.2426036 -2.3912496e-15]
 [   0.076923077      2.3136095      4.5502959      4.5976331      4.2071006      4.0769231      4.0769231      4.5384615      5.4615385      5.9230769      5.4615385      4.6863905      2.3431953 -4.5092135e-15]
 [    0.61538462      2.3550296      4.0946746      4.4733728      4.5798817      4.6153846      4.6153846      5.0769231              6      6.4615385              6       5.183432       2.591716 -4.9874634e-15]
 [     1.3076923      2.4792899      3.6508876      4.2307692      4.6923077      4.8934911       4.964497      5.3550296      6.0650888       6.408284      6.0177515      5.2544379      2.6272189 -5.0557849e-15]
 [     2.3846154      2.8106509      3.2366864      3.6923077      4.1538462      4.5207101      4.8402367      4.9822485      4.9467456      4.8757396      4.7337278       4.260355      2.1301775  -4.099285e-15]
 [     2.5384615      2.6449704      2.7514793      3.0650888      3.4201183      3.8461538      4.3076923      4.3076923      3.8461538      3.5029586      3.3964497      3.0532544      1.5266272 -2.9378209e-15]
 [     1.4615385       1.816568      2.1715976      2.3195266      2.4260355      2.7692308      3.2307692      3.2307692      2.7692308      2.3431953      1.9881657      1.5621302     0.78106509 -1.5030712e-15]
 [    0.69230769      1.1538462      1.6153846      1.5739645      1.4319527      1.6923077      2.1538462      2.2248521      1.9053254      1.5384615      1.0769231     0.63905325     0.31952663 -6.1489275e-16]
 [    0.15384615     0.61538462      1.0769231     0.82840237     0.43786982     0.61538462      1.0769231      1.2721893      1.2011834              1     0.53846154     0.14201183    0.071005917 -1.3664283e-16]
 [             0     0.28402367     0.56804734     0.37869822    0.094674556     0.18934911     0.47337278     0.61538462     0.61538462     0.52071006     0.23668639              0              0              0]
 [             0    0.035502959    0.071005917    0.047337278     0.01183432    0.023668639    0.059171598    0.076923077    0.076923077    0.065088757    0.029585799              0              0              0]
 [             0              0              0              0              0              0              0              0              0              0              0              0              0              0]
 [             0              0              0              0              0              0              0              0              0              0              0              0              0              0]
]

lis · 2008-05-28 17:25 @767

Hola, muchas gracias por responderme.

Yo descargo datos grib, los que transformo a .ctl. Estos datos tienen 27 variables (temperatura, viento, presion, etc.), entonces para trabajar con una de ellas uso el wgrib y en este caso ocupo la temperatura como ejemplo. En Perl hago lo siguiente:

Syntax: [ Download ] [ Hide ]

Using perl Syntax Highlighting

#!/usr/bin/perl

open (TAREA, ">tarea6.dat");

$file_save="tarea6.dat";

$filename = "/home/lis/grads/tarea6/cl_wafsgfs_P_t18z_intdsk12";

open(GUARDAR,"$file_save");

system"(wgrib $filename |grep 'kpds5=11' |grep 'kpds6=105' |grep 'kpds7=2' |wgrib -i $filename -text -o $file_save)";

close(GUARDAR);

Y esto me genera un archivo llamado tarea.dat que contiene una columna de 2450 datos, que en este caso son de temperatura. Lo que tengo que hacer es ingresar una coordenada (latitud, longitud) y me tiene que mostrar los 4 números que rodean a esta coordenada y luego ingresarlos a la ecuación de interpolación lineal y que me muestre resultado de la temperatura interpolada.

Gracias por ayudarme, y estoy muy contenta por la pronta respuesta.

explorer · 2008-05-31 19:37 @859

¿2450 valores? Pero si a mi me salen 3447...

¿De qué directorio ftp del NOAA te bajas los ficheros?

lis · 2008-06-04 01:16 @095

Hola explorer. Gracias por tu ayuda. Me ha sido de gran utilidad.

Mis datos eran 2450 porque yo acoto los rangos a mi medida, es decir, mi grilla es: -70° latitud sur a -8.75° latitud sur y -120° longitud oeste a -60° longitud oeste.

Gracias a ti ya resolví esa parte, ahora estoy haciendo un script en Perl con la interpolacióon de punto más cercano o vecino más próximo (nearest neighbor), y mi duda es la siguiente:

Con la interpolación bilineal (hecha en el script anterior) yo ingresaba dos datos que por ejemplo son: (x.y) = (longitud, latitud) = (-70.73, 30.25) y me muestra los 4 puntos que rodean al punto ingresado (x,y), con sus respectivas coordenadas y como yo estaba ocupando la temperatura, también me da la temperatura de esos 4 puntos y por último me da la temperatura interpolada bilinealmente.

Ahora para esa misma coordenada ingresada (-70.73, -30.25), en Perl tengo que aplicarle la interpolación del punto más cercano, y como resultado me tiene que mostrar las coordenadas del punto más cercano a (-70.73, -30.25) y la temperatura ahora interpolada por la interpolación del punto más cercano.

Recurro a ti nuevamente porque para esta nueva interpolación no hay una ecuación, y quiero saber si hay algún comando en Perl que haga esto de encontrar el punto más cercano de cualquier variable.

Cualquier ayuda te lo agradecería mucho, Explorer.

explorer · 2008-06-05 06:30 @312

Bien, tu rejilla es más pequeña, pero habrás tenido en cuenta que lo que se guarda en un grib no es una rejilla cartesiana, sino que en cada latitud se guarda un número diferente de muestras (menos cuanto más te acercas al polo). Hay que tener cuidado con eso cuando se extraen los datos.

En cuanto a lo del comando Perl, no existe algo parecido. Hay que realizar las cuentas a mano. Ya lo has visto en otro hilo.

Un detalle importante. No es lo mismo hacer esto para uno o dos puntos, que para millones de puntos. Perl no es el lenguaje más adecuado para realizar grandes procesamientos de cálculo, ya que entonces es mejor usar soluciones compiladas. En ese caso se hace imprescindible usar una librería matemática, como la PDL.

Publicidad

Socializa