/* * Copyright (C) 1996-2011 by ECOLE DES MINES DE PARIS * * Permission is hereby granted, free of charge, to any person obtaining a copy * of this software and associated documentation files (the "Software"), to deal * in the Software without restriction, including without limitation the rights * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell * copies of the Software, and to permit persons to whom the Software is * furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN * THE SOFTWARE. */ /*-------------------------------------------------------------------------*/ /* ECOLE DES MINES DE PARIS */ /* CENTRE D'ENERGETIQUE - GROUPE TELEDETECTION & MODELISATION */ /* Rue Claude Daunesse, BP 207 */ /* 06904 Sophia Antipolis cedex, FRANCE */ /* Tel (+33) 04 93 95 74 49 Fax (+33) 04 93 95 75 35 */ /* E-mail : (name)@cenerg.cma.fr */ /*-------------------------------------------------------------------------*/ /* meteosatlib.c */ /* December 1996. L. Wald & M. Lefevre */ /* Some useful functions for the geometry of the satellite, and */ /* conversion of latitude-longitude into image coordinates and vice-versa. */ /*-------------------------------------------------------------------------*/ #include int latlon_to_HR(); int latlon_to_B2(); int HR_to_latlon(); int B2_to_latlonHRb2(); int B2_to_latlonHRb0(); int HR_to_HRb2(); int HR_to_B2(); int HR_to_HRb2v(); int HRb2_to_B2(); int B2_to_HRb2(); int satelev(); int satazim(); /****************************************************************************** Meteosat images treated in this library are 2 products : 1. a high resolution image (above noted HR) with 2500 lines and 2500 pixels per line oriented NW (lig 1, col 1) to SE (lig 2500, col 2500). The pixel size is 5 km. When such coordinates in lines and columns are used, they always refer to the image correctly oriented NW-SE. The acquired image is actually inverse, taken from SE to NW, so that the first pixel acquired is the south-easternmost pixel on the image (2500, 2500). 2. a B2 image (above noted B2) with 416 lines and 416 pixels per line, subsampled from the high resolution as follow : starting from the first pixel acquired (2500, 2500) are taken each pixel every 6 lines and every 6 columns. These selected pixels are noted HRb2 in the high resolution product, and B2 in the B2 product. Thus a pixel B2 represents a pack of 6x6 HR pixels whose down- right one substitutes the other ones. It is noted that 2500 cannot be divided by 6, so that the last 4 lines and columns in the acquired direction are ignored by the B2 product. The first HRb2 pixel is positioned (10,10) on the high resolution image, while the last one is (2500, 2500). This library provides functions to calculate position of pixels from one product to the other, as well as conversion between pixel positions and geographical coordinates (lat, lon). *******************************************************************************/ /* ***** M E T E O S A T P I X E L C O N V E R S I O N ***** */ /* ***** B 2 / H I G H _ R E S O L U T I O N ***** */ /******************************************************************************/ int HR_to_HRb2(ligHR, colHR, ligHRb2, colHRb2) int ligHR, colHR; int *ligHRb2, *colHRb2; /* The procedure 'HR_to_HRb2' provides for any pixel HR (lig, col in the NW-SE oriented high resolution image, the corresponding pixel HRb2 (lig, col in NW-SE oriented high resolution image) used for the B2 format. The return value is : 1 in case of wrong input data (i.e. pixel outside the HR image frame), 0 if ok, and when different from 0, lig and col are set to -999. */ { *ligHRb2 = -999 ; *colHRb2 = -999 ; if(ligHR < 1 || ligHR > 2500) return 1; if(colHR < 1 || colHR > 2500) return 1; if(ligHR < 5 || colHR < 5) return 2; *ligHRb2 = ligHR; *colHRb2 = colHR; while(((*ligHRb2 - 4) % 6) != 0) *ligHRb2 = *ligHRb2 + 1 ; while(((*colHRb2 - 4) % 6) != 0) *colHRb2 = *colHRb2 + 1 ; return 0; } int HRb2_to_B2(ligHRb2, colHRb2, ligB2, colB2) int ligHRb2, colHRb2; int *ligB2, *colB2; /* The procedure 'HRb2_to_B2' provides for any pixel HRb2 (lig, col in the NW-SE oriented high resolution image), the corresponding pixel in the B2 format (lig, col in NW-SE oriented B2 image). The return value is : 1 in case of wrong input data (i.e. pixel outside the HR image frame), 4 in case of wrong input data (i.e. pixel not used for B2 image), 0 if ok, and when different from 0, lig and col are set to -999. */ { *ligB2 = -999 ; *colB2 = -999 ; if(ligHRb2 < 1 || ligHRb2 > 2500) return 1; if(colHRb2 < 1 || colHRb2 > 2500) return 1; if(ligHRb2 < 10 || colHRb2 < 10) return 4; if( ((ligHRb2 - 4) % 6 ) * ((colHRb2 - 4) % 6 ) != 0 ) return 4 ; *ligB2 = (ligHRb2 - 4) / 6 ; *colB2 = (colHRb2 - 4) / 6 ; return 0; } int B2_to_HRb2(ligB2, colB2, ligHRb2, colHRb2) int ligB2, colB2; int *ligHRb2, *colHRb2; /* The procedure 'B2_to_HRb2' provides for any pixel B2 (lig, col in NW-SE oriented B2 image), the corresponding pixel HRb2 (lig, col in the NW-SE oriented high resolution image. The return value is : 1 in case of wrong input data (i.e. pixel outside the B2 image frame), 0 if ok, and when different from 0, lig and col are set to -999. */ { *ligHRb2 = -999 ; *colHRb2 = -999 ; if(ligB2 < 1 || ligB2 > 416) return 1; if(colB2 < 1 || colB2 > 416) return 1; *ligHRb2 = 4 + ligB2 * 6; *colHRb2 = 4 + colB2 * 6; return 0; } int HR_to_B2(ligHR, colHR, ligB2, colB2) int ligHR, colHR; int *ligB2, *colB2; /* The procedure 'HR_to_B2' provides for any pixel HR (lig, col in NW-SE oriented high resolution image), the corresponding pixel in the B2 format (lig, col in NW-SE oriented B2 image). The return value is : 1 in case of wrong input data (i.e. pixel outside the HR image frame), 0 if ok, and when different from 0, lig and col are set to -999. */ { int ligHRb2, colHRb2; int i; *ligB2 = -999 ; *colB2 = -999 ; i = HR_to_HRb2(ligHR, colHR, &ligHRb2, &colHRb2); if(i != 0) return i; i = HRb2_to_B2(ligHRb2, colHRb2, ligB2, colB2); return i; } /******************************************************************************/ /* ***** G E O _ L O C A T I O N O F M E T E O S A T P I X E L S ***** */ /* ***** B 2 & H I G H _ R E S O L U T I O N ***** */ /******************************************************************************/ int latlon_to_HRb2v(lat, lon, ligHRb2v, colHRb2v, usable, dist) double lat, lon; int ligHRb2v[], colHRb2v[]; int usable[]; float dist[]; /* The procedure 'latlon_to_HRb2v' provides for any geographical coordinates (degrees, positive NE and negative SW) the sixteen nearest pixels HRb2 (lig, col in NW-SE oriented high resolution image) used for the B2 image, and sort by increasing distances with the HR pixel corresponding to the geographical coordinates. The parameter usable is 1 if the resulting pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The return value is : 1 in case of wrong input geographical coordinates, 2 in case of no B2 corresponding pixels, 3 in case of non usable pixel, 0 if OK, and when different from 0 lig, col and distances are set to -999 */ { int ligHR, colHR; int ligHRb2, colHRb2; int i, j, k, a2, b2; float y; for( i=0 ; i<16 ; i=i+1 ) { ligHRb2v[i] = -999; colHRb2v[i] = -999; dist[i] = -999.0; usable[i] = 0; } i = latlon_to_HR(lat, lon, &ligHR, &colHR, usable); if(i != 0) return i; i = HR_to_HRb2(ligHR, colHR, &ligHRb2, &colHRb2); if(i != 0) return i; for( i=-2, k=0 ; i<=1 ; i=i+1 ) for( j=-2 ; j<=1 ; j=j+1, k=k+1 ) { ligHRb2v[k] = ligHRb2 + (i * 6); colHRb2v[k] = colHRb2 + (j * 6); usable[k] = 1; if(ligHRb2v[k] < 10 || colHRb2v[k] < 10) usable[k] = 0; if(ligHRb2v[k] > 2500 || colHRb2v[k] > 2500) usable[k] = 0; if(usable[k] == 0) continue; a2 = (ligHRb2v[k] - ligHR) * (ligHRb2v[k] - ligHR); b2 = (colHRb2v[k] - colHR) * (colHRb2v[k] - colHR); dist[k] = sqrt( a2 + b2 ); } for( i=0 ; i<16 ; i=i+1 ) { for( j=i+1 ; j<16 ; j=j+1 ) if(dist[j] < dist[i]) { y = dist[i]; dist[i] = dist[j]; dist[j] = y; k = ligHRb2v[i]; ligHRb2v[i] = ligHRb2v[j]; ligHRb2v[j] = k; k = colHRb2v[i]; colHRb2v[i] = colHRb2v[j]; colHRb2v[j] = k; k = usable[i]; usable[i] = usable[j]; usable[j] = k; } /* printf("v%-2d. %4d x %4d d =%6.1f usable = %d\n",i,ligHRb2v[i],colHRb2v[i],dist[i], usable[i]); */ } return 0; } int latlon_to_HR(lat,lon,ligHR,colHR,usable) double lat,lon ; int *ligHR,*colHR,*usable ; /* The procedure 'latlon_to_HR' provides the geographical coordinates (degrees, positive NE and negative SW) for any pixel HR (lig, col in NW-SE oriented high resolution image). The parameter usable is 1 if the resulting pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The returned value is : 1 in case of wrong input geographical coordinates, 3 in case of non usable pixel, 0 if OK, and when different from 0 lig and col are set to -999. */ { double tlat,tlon,tlon0 ; double sinlat,coslat,tanlat,sinlon,coslon,sinlon0,coslon0; double ReRp,RpdRe,Hs2,Rs,Rp2,Re2; double teta,rom,rom2,r1,y2 ; double vt,wt,xt,yt,zt,costeta,sinteta,px,py ; double xr,yr ; const double deg_rad = 0.0174533 ; /* converts decimal degrees into radians */ const double Re = 6378.15 ; /* earth equatorial radius */ const double Rp = 6356.75 ; /* earth polar radius */ const double Hs = 35786.0 ; /* altitude of satellite */ const double lon0 = 0.0 ; /* longitude of satellite (Meteosat) */ const double deltax = 0.0072 ; /* radiometer scanning steps (for Meteosat) : */ const double deltay = 0.0072 ; /* E-W = 18.0/2500 (col) N-S = 18.0/2500 (lig) */ *ligHR = -999; *colHR = -999; *usable = 0; if(lat < -90.0 || lon < -90) return 1; if(lat > 90.0 || lon > 90) return 1; tlat = lat * deg_rad ; tlon = lon * deg_rad ; coslat = cos(tlat) ; sinlat = sin(tlat) ; coslon = cos(tlon) ; ReRp = Re * Rp ; Hs2 = Hs * Hs ; Rp2 = Rp * Rp ; Re2 = Re * Re ; rom = ReRp / sqrt(Rp2 * coslat * coslat + Re2 * sinlat * sinlat) ; rom2 = rom * rom ; y2 = Hs2 + rom2 - 2 * Hs * rom * coslat * coslon ; r1 = y2 + rom2 ; if (r1 > Hs2) return 3; RpdRe = Rp / Re ; Rs = Re + Hs ; tanlat = tan(tlat) ; sinlon = sin(tlon) ; tlon0 = lon0 * deg_rad ; sinlon0 = sin(tlon0) ; coslon0 = cos(tlon0) ; teta = atan(RpdRe * tanlat) ; costeta = cos(teta) ; sinteta = sin(teta) ; xt = Re * costeta * coslon ; yt = Re * costeta * sinlon ; zt = Rp * sinteta ; wt = xt - Rs * coslon0 ; vt = yt - Rs * sinlon0 ; px = atan(( coslon0 * vt - wt * sinlon0 ) / ( sinlon0 * vt + wt * coslon0) ) ; py = atan(( zt * (tan(px) * sinlon0 - coslon0) / wt ) * cos(px)) ; px = px / deg_rad ; py = py / deg_rad ; xr = px / deltax ; yr = py / deltay ; /* MODIF 15/06/2009 pour corriger la non reversibilite latlon_to_ligcol if (xr >= 0.0) xr = floor(px/deltax)+0.5 ; if (xr < 0.0) xr = floor(px/deltax)-0.5 ; if (yr >= 0.0) yr = floor(py/deltay)+0.5 ; if (yr < 0.0) yr = floor(py/deltay)-0.5 ; */ xr = floor(px/deltax)+0.5 ; yr = floor(py/deltax)+0.5 ; *colHR = (int) floor(xr+1250.5); *ligHR = (int) floor(yr+1250.5); *ligHR = 2500 - *ligHR +1; *colHR = 2500 - *colHR +1; *usable = 1 ; return 0 ; } int HR_to_latlon(ligHR, colHR, lat, lon, usable) int ligHR, colHR, *usable ; double *lat, *lon ; /* The procedure 'HR_to_latlon' provides for any pixel HR (lig, col in NW-SE oriented high resolution image) its geographical coordinates (degrees, positive NE and negative SW). The parameter usable is 1 if the input pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The returned value is : 1 in case of wrong input data (i.e. lig and col outside the image frame), 3 in case of non usable pixel, 0 if OK, and when different from 0 lat and lon are set to -999.0 */ { const double deg_rad = 0.0174533 ; /* converts decimal degrees into radians */ const double Re = 6378.15 ; /* earth equatorial radius */ const double Rp = 6356.75 ; /* earth polar radius */ const double Hs = 35786.0 ; /* altitude of satellite */ const double Rs = Re + Hs ; /* orbit radius */ const double yk2 = Rs * Rs / (Re * Re) ; const double A = 1.0/297.0; /* earth oblatness */ const double lon0 = 0.0 ; /* longitude of satellite (Meteosat) */ const double deltaX = 0.0072 ; /* Meteosat radiometer scanning steps : */ const double deltaY = 0.0072 ; /* E-W = 18.0/2500 (col) N-S = 18.0/2500 (lig) */ static double X,Y,cosX,tanX,tan2X,tanY,tan2Y; /* coordinates of point P(X,Y) */ static double tlon0, coslon0,sinlon0; /* coordinates of Meteosat */ static double val1,val2,vmu,xt,yt,zt,teta; *lat = -999.0; *lon = -999.0; *usable = 0; if(ligHR < 1 || ligHR > 2500) return 1; if(colHR < 1 || colHR > 2500) return 1; X = 2500 - colHR + 1; /* returns image from NW-SE to SE-NW acquisition */ Y = 2500 - ligHR + 1; X = X - 1250.5; /* translates subsatellite point (1250,1250) to (0,0) */ Y = Y - 1250.5; X = X * deltaX * deg_rad; Y = Y * deltaY * deg_rad; cosX = cos(X); tanX = tan(X); tan2X = tanX*tanX; tanY = tan(Y); tan2Y = tanY*tanY; tlon0 = lon0 * deg_rad; coslon0 = cos(tlon0); sinlon0 = sin(tlon0); val1=1.0 + tan2X; val2=1.0 + tan2Y*(A+1)*(A+1); if( (val1 * val2) > (yk2 / (yk2 - 1)) ) return 3; vmu = ( Rs - Re * sqrt(yk2 - (yk2 - 1) * val1 * val2) ) / (val1 * val2); xt = (Rs * coslon0) + ( vmu * (tanX * sinlon0 - coslon0) ); yt = (Rs * sinlon0) - ( vmu * (tanX * coslon0 + sinlon0) ); zt = vmu * tanY / cosX; teta=asin(zt / Rp); *lat = atan( tan(teta) * Re / Rp) / deg_rad; *lon = atan(yt/xt) / deg_rad; *usable = 1; return 0; } int B2_to_latlonHRb2(ligB2, colB2, lat, lon, usable) int ligB2, colB2, *usable ; double *lat, *lon ; /* The procedure 'B2_to_latlonHRb2' provides for any pixel B2 (lig, col in NW-SE oriented B2 image) the geographical coordinates (degrees, positive NE and negative SW) of the corresponding HRb2 pixel whose value in the high resolution image is used for the B2 format. The parameter usable is 1 if the input pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The returned value is : 1 in case of wrong input data (i.e. lig col outside the B2 image frame), 3 in case of non usable pixel, 0 if OK, and when different from 0 lat and lon are set to -999.0 */ { int i; int ligHRb2, colHRb2; *lat = -999.0; *lon = -999.0; *usable = 0; i = B2_to_HRb2(ligB2, colB2, &ligHRb2, &colHRb2); if(i != 0 ) return i; i = HR_to_latlon(ligHRb2, colHRb2, lat, lon, usable); return i; } int B2_to_latlonHRb0(ligB2, colB2, lat, lon, usable) int ligB2, colB2, *usable ; double *lat, *lon ; /* The procedure 'B2_to_latlonHRb0' provides for any pixel B2 (lig, col in NW-SE oriented B2 image).the geographical coordinates (degrees, positive NE and negative SW) of the corresponding HRb0 pixel in the high resolution image. HRb0 is the virtual pixel centered in the 6x6 pixels whose values are substituted by the right down value. Its coordinates are interpolated from two surrounding pixels. The parameter usable is 1 if the input pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The returned value is : 1 in case of wrong input data (i.e. lig col outside the B2 image frame), 3 in case of non usable pixel, 0 if OK, and when different from 0 lat and lon are set to -999.0 */ { int i; int ligHRb2, colHRb2; int ligHRb0, colHRb0; int ligHRb00, colHRb00; double lat0, lon0; double lat00, lon00; *lat = -999.0; *lon = -999.0; *usable = 0; i = B2_to_HRb2(ligB2, colB2, &ligHRb2, &colHRb2); if(i != 0) return i; ligHRb0 = ligHRb2 - 3; colHRb0 = colHRb2 - 3; ligHRb00 = ligHRb2 - 2; colHRb00 = colHRb2 - 2; i = HR_to_latlon(ligHRb0, colHRb0, &lat0, &lon0, usable); if(i != 0) return i; i = HR_to_latlon(ligHRb00, colHRb00, &lat00, &lon00, usable); if(i != 0) return i; *lat = (lat0 + lat00) / 2.0 ; *lon = (lon0 + lon00) / 2.0 ; return 0; } int latlon_to_B2(lat, lon, ligB2, colB2, usable) double lat, lon ; int *ligB2, *colB2, *usable ; /* The procedure 'latlon_to_B2' provides for any geographical coordinates (degrees, positive NE and negative SW), the corresponding pixel B2 (lig, col in NW-SE oriented B2 image). The parameter usable is 1 if the pixel is located in the useful part of the Meteosat image, 0 otherwise (i.e. not in the field of view). The returned value is : 1 in case of wrong input geographical coordinates, 3 in case of non usable pixel, 0 if OK, and when different from 0 lig and col are set to -999 */ { int i; int ligHR, colHR; *ligB2 = -999; *colB2 = -999; *usable=0; i=latlon_to_HR(lat, lon, &ligHR, &colHR, usable); if(i != 0) return i; i=HR_to_B2(ligHR, colHR, ligB2, colB2); return i; } /******************************************************************************/ /* **** G E O M E T R I C P A R A M E T E R S ***** */ /******************************************************************************/ int satelev(lat, lon, lon0, gamma) double lat, lon, lon0 ; double *gamma ; /* The procedure 'satelev' provides the elevation of the satellite above */ /* horizon. Returns 0 if OK, 1 otherwise. lat is the latitude, in decimal */ /* degree, positive north. lon is the longitude, in decimal degree, posi- */ /* tive east. lon0 is the longitude of the satellite, same conventions than */ /* lon. gamma is the elevation, in radian. */ { double tlat, tlon, tlon0, sin_gamma ; double x, y, z; const double deg_rad = 0.0174533 ; /* converts decimal degrees into radians*/ const double Re = 6378.15 ; /* earth equatorial radius */ const double Hs = 35786.0 ; /* altitude of satellite */ tlat = lat * deg_rad ; tlon = lon * deg_rad ; tlon0 = lon0 * deg_rad ; x = cos(tlat) ; y = cos(tlon - tlon0) ; z = Re + Hs ; sin_gamma = ( z * x * y - Re ) / sqrt( Hs * Hs + 2 * Re * z * (1 - x * y) ); *gamma = asin(sin_gamma); return 0 ; } int satazim(lat, lon, phiv) double lat, lon; double *phiv; /* The procedure 'satazim' provides the satellite azimuth angle */ /* Returns 0 if OK, 1 otherwise. lat is the latitude, in decimal degree */ /* positive north. lon is the longitude, in decimal degree, positive east */ /* lon0 is the longitude of the satellite, same conventions than */ /* lon. phiv is the azimuthal angle, in radian. */ { double tlat, tlon; const double pi = 3.141593 ; const double deg_rad = 0.0174533 ; /* converts decimal degrees into radians*/ tlat = lat * deg_rad ; tlon = lon * deg_rad ; if(sin(tlat) == 0) *phiv = pi; else *phiv = (atan( tan(tlon) / sin(tlat)) + pi); return 0; }