/* * Copyright (C) 2004-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 */ /*--------------------------------------------------------------------------------------------------------------------*/ /* HELIOCLIM : msglib.c Feb.2004 L. Wald & M. Lefevre HELIOCLIM3 : msglib.v2.c 2008 Some useful functions for the geometry of the satellite MSG : Conversion of geographic coordinates into image coordinates and vice-versa. For all images lines and columns start from 1 and incrementation is NE-SW oriented. Conversion of high resolution (3km, 3712 x 3712 pixels)) image coordinates into low resolution (10km, 1238 x 1238 pixels) and vice-versa Note that the radiometer has a total scanning grid of 3750 x 3750 lig and col with a scan step of 251.53 microrad (=0.0144116°) for 3 lines (sources: EUMETSAT MSG Level 1.5 Image Data Format Description, Issue 2, nov 2001,p.169) It follows that 3750 lines are scanned within 18.0145° field of view angle. M S G P I X E L C O N V E R S I O N & G E O L O C A T I O N ----------------------------------------------------------------------------------------------------------------------*/ int latlon_to_ligcol(); int ligcol_to_latlon(); int I3_to_latlon(); int latlon_to_I3(); int I10_to_latlon(); /* supprimer ? */ int latlon_to_I10(); /* supprimer ? */ int I10_to_I3(); /* supprimer ? */ int I3_to_I10(); /* supprimer ? */ int satelev(); int satazim(); int hc2delta_t(); double hc3delta_t(); void load_MSG_info(); /* AJOUT 2008 msglib.v2.c lecture des parametres du satellite */ extern void end(); /*---------------------------------------------------------------------------------------------------------------------- Provides the geographical coordinates (degrees, positive NE and negative SW) for any pixel (lig, col) from (1,1) to (I3_Nblig,I3_Nbcol) in NW-SE oriented high resolution image). The parameter usable is 1 if the resulting pixel is located in the useful part of the 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 latlon_to_ligcol(deltaX,deltaY,nblig,nbcol,lat,lon,lig,col,usable) double deltaX,deltaY; int nblig,nbcol; double lat,lon ; int *lig,*col,*usable ; { 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 ; *lig = -999; *col = -999; *usable = 0; if(lat < -90.0 || lon < -90) return 1; if(lat > 90.0 || lon > 90) return 1; tlat = lat * deg2rad ; tlon = lon * deg2rad ; 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 * deg2rad ; 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 / deg2rad ; py = py / deg2rad ; xr = px / deltaX ;/* coordinates centered over subsatellite point */ yr = py / deltaY ; /* MODIF 28/05/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 ; si lon>0 alors col decalee de +1 if (yr >= 0.0) yr = floor(py/deltaY)+0.5 ; if (yr < 0.0) yr = floor(py/deltaY)-0.5 ; si lat<0 alors lig decalee de +1 */ xr = floor(px/deltaX)+0.5 ; yr = floor(py/deltaY)+0.5 ; *col = (int) floor(xr+(nbcol/2)+0.5); *lig = (int) floor(yr+(nblig/2)+0.5); *lig = nblig - *lig +1 ; /* renversement N-S de l'image */ *col = nbcol - *col +1 ; *usable = 1 ; return 0 ; } /*---------------------------------------------------------------------------------------------------------------------- Provides for any pixel (lig, col) from (1,1) 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 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 */ int ligcol_to_latlon(deltaX,deltaY,nblig,nbcol,lig,col,lat,lon,usable) double deltaX,deltaY; int nblig,nbcol,lig,col,*usable ; double *lat,*lon ; { const double Rs = Re + Hs ; /* orbit radius */ const double yk2 = Rs * Rs / (Re * Re) ; const double A = 1.0/297.0; /* earth oblatness */ double tlon0, coslon0,sinlon0; /* coordinates of MSG */ double X,Y,cosX,tanX,tan2X,tanY,tan2Y;/* coordinates of point(X,Y)*/ double val1,val2,vmu,xt,yt,zt,teta; *lat = -999.0; *lon = -999.0; *usable = 0; if(lig < 1 || lig > nblig) return 1; if(col < 1 || col > nbcol) return 1; X = nbcol - col + 1; /* returns image from NW-SE to SE-NW acquisition*/ Y = nblig - lig + 1; X = X - (nbcol/2) - 0.5;/* translates subsatellite point to (0,0) */ Y = Y - (nblig/2) - 0.5; X = X * deltaX * deg2rad; Y = Y * deltaY * deg2rad; cosX = cos(X); tanX = tan(X); tan2X = tanX*tanX; tanY = tan(Y); tan2Y = tanY*tanY; tlon0 = Lon0 * deg2rad; /* longitude du satellite */ coslon0 = cos(tlon0); sinlon0 = sin(tlon0); val1=1.0 + tan2X; val2=1.0 + tan2Y*(A+1)*(A+1); vmu = yk2 - (yk2 - 1)*val1*val2; if( vmu < 0 ) return 3; /* exclude a negative argument for sqrt() */ vmu = ( Rs - Re * sqrt(vmu) ) / (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) * rad2deg; *lon = atan(yt/xt) * rad2deg; *usable = 1; return 0; } /*--------------------------------------------------------------------------------------------------------------------*/ int I3_to_latlon(lig, col, lat, lon, usable) int lig, col, *usable ; double *lat, *lon ; { static double deltaX = Scanstep/3.0*rad2deg; /* 1 pixel scanned / image col */ static double deltaY = Scanstep/3.0*rad2deg; /* ............... / lig */ return ligcol_to_latlon(deltaX,deltaY,I3_NBLIG,I3_NBCOL,lig,col,lat,lon,usable); } /*--------------------------------------------------------------------------------------------------------------------*/ int latlon_to_I3(lat, lon, lig, col, usable) double lat,lon ; int *lig,*col,*usable ; { static double deltaX = Scanstep/3.0*rad2deg; /* 1 pixel scanned / image col */ static double deltaY = Scanstep/3.0*rad2deg; /* ............... / lig */ return latlon_to_ligcol(deltaX,deltaY,I3_NBLIG,I3_NBCOL,lat,lon,lig,col,usable); } /*--------------------------------------------------------------------------------------------------------------------*/ int I10_to_latlon(lig, col, lat, lon, usable) int lig, col, *usable ; double *lat, *lon ; { int ier=0; const double deltaX = Scanstep*rad2deg; /* 3 pixels scanned /col */ const double deltaY = Scanstep*rad2deg; /* ................ /line */ ier=ligcol_to_latlon(deltaX,deltaY,I10_NBLIG,I10_NBCOL,lig,col,lat,lon,usable); return ier; } /*--------------------------------------------------------------------------------------------------------------------*/ int latlon_to_I10(lat, lon, lig, col, usable) double lat,lon ; int *lig,*col,*usable ; { int ier=0; const double deltaX = Scanstep*rad2deg; /* 3 pixels scanned /col */ const double deltaY = Scanstep*rad2deg; /* ................ /line */ ier=latlon_to_ligcol(deltaX,deltaY,I10_NBLIG,I10_NBCOL,lat,lon,lig,col,usable); return ier; } /*---------------------------------------------------------------------------------------------------------------------- Provides for any pixel in low resolution the corresponding pixel in the high resolution image (top left of 3x3 pixels). The return value is : 1 in case of wrong input data (i.e. pixel outside the image frame), 0 if ok, and when different from 0, lig and col are set to -999. */ int I10_to_I3(lig10, col10, lig3, col3) int lig10, col10; int *lig3, *col3; { *lig3 = -999 ; *col3 = -999 ; if(lig10 < 1 || lig10 > I10_NBLIG) return 1; if(col10 < 1 || col10 > I10_NBCOL) return 1; *lig3 = lig10 * 3 - 2; *col3 = col10 * 3 - 2; return 0; } /*---------------------------------------------------------------------------------------------------------------------- Provides for some pixels (top left of 3x3 pixels in the high resolution) the corresponding pixel in the low resolution image. The return value is : 1 in case of wrong input data (i.e. pixel outside the image frame), 0 if ok, and when different from 0, lig and col are set to -999. */ int I3_to_I10(lig3, col3, lig10, col10) int lig3, col3; int *lig10, *col10; { *lig10 = -999 ; *col10 = -999 ; if(lig3 < 1 || lig3 > I3_NBLIG) return 1; if(col3 < 1 || col3 > I3_NBCOL) return 1; *lig10=lig3; *col10=col3; /* lig and col should be centered on 3x3 pixels ??? revoir */ if((*lig10 % 3 )==1) *lig10=*lig10+1; else if(*lig10 % 3 == 0) *lig10=*lig10-1; if((*col10 % 3 )==1) *col10=*col10+1; else if(*col10 % 3 == 0) *col10=*col10-1; *lig10 = 1 + *lig10 / 3 ; *col10 = 1 + *col10 / 3 ; return 0; } /*---------------------------------------------------------------------------------------------------------------------- Provides the elevation gamma (rd) of the satellite above horizon. */ int satelev(lat, lon, gamma) double lat, lon; /* decimal degrees positive N and E */ double *gamma; /* elevation in radian */ { double tlat, tlon, tlon0, sin_gamma, x, y, z; tlat = lat * deg2rad ; tlon = lon * deg2rad ; tlon0 = Lon0 * deg2rad ;/* longitude of satellite */ 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 ; } /*---------------------------------------------------------------------------------------------------------------------- Provides the satellite azimuth angle phiv (rd). */ int satazim(lat, lon, phiv) double lat, lon; /* decimal degrees positive N and E */ double *phiv; /* azimuthal angle, in radian */ { double tlat, tlon; tlat = lat * deg2rad ; tlon = lon * deg2rad ; if(sin(tlat) == 0) *phiv = PI; else *phiv = (atan( tan(tlon) / sin(tlat)) + PI); return 0; } /*---------------------------------------------------------------------------------------------------------------------- Décalage entre l'heure de début de balayage de la ligne (résolution 10 km) et l'heure de début d'acquisition de l'image haute resolution = 12.5'/3750/3 lig balayées 3x3 Returns 0 if OK, 1 otherwise. */ int hc2delta_t(lig10, delta_t) int lig10; /* lines counted from 1 in the low resolution */ double *delta_t; /* in decimal hour */ { static const int offlig=19; static double nbstep=0; if(nbstep==0) nbstep=(double)(I3_NBLIG + 2 * offlig) / 3.0; *delta_t = -999.0; if(lig10<=0 || lig10>I10_NBLIG) return 1; *delta_t = (double)(I10_NBLIG-lig10+1+(offlig/3.0));/* inversion N-S */ *delta_t*= 12.5 / (double)nbstep / 60.0; return 0; } /*---------------------------------------------------------------------------------------------------------------------- Décalage entre l'heure de début de balayage de la ligne et l'heure de début d'acquisition de l'image haute resolution = 12.5'/3750/3 lig balayées 3 par 3 Return delta_t in decimal hour */ double hc3delta_t(lig) int lig; /* lines counted from 1 in geogr. grid (N-S) in the high resolution */ { static const int offlig=19; static double nbstep=0; double delta_t=-999.0; /* in decimal hour */ if(nbstep==0) nbstep=(double)(I3_NBLIG + 2 * offlig) / 3.0; if(lig>0 && lig<=I3_NBLIG) { delta_t = (double)(I3_NBLIG-lig+1+offlig)/3.0;/* inversion N-S */ if(0) printf("# delta_t=%5.2f'\n",delta_t*12.5/(double)nbstep); delta_t*= 12.5 / (double)nbstep / 60.0; } return delta_t; } /*---------------------------------------------------------------------------------------------------------------------- Chargement des parametres du satellite */ void load_MSG_info(yy0,mm0,dd0,hr0,mn0,mysat,a1,b1,a2,b2) int yy0,mm0,dd0,hr0,mn0,*mysat; /* date et satellite operationnel */ double *a1,*b1,*a2,*b2; /* parametres de calibration des canaux 1 et 2 */ { int yyyymmdd,hr,mn; char *buf=(char *) malloc(200); FILE *fi=fopen(INFO_MSG,"r"); int date= (yy0 + (yy0 <50 ? 2000:0))*10000+mm0*100+dd0; *mysat=0; *a1=*a2=*b1=*b2=0.0; if(yy0<0 || mm0<0 || dd0<0 || hr0<0 || mn0<0 || date < 20040201) return; while(!feof(fi)){ memset(buf,'#',strlen(buf)); if(fgets(buf,200,fi)==NULL) break; if(*buf=='#' || *buf=='\n') continue; if(sscanf(buf,"%d %02d:%02d %d %lf %lf %lf %lf",&yyyymmdd,&hr,&mn,mysat,a1,b1,a2,b2)!=8) end(2,INFO_MSG); if(date > yyyymmdd) break; if(date == yyyymmdd) { if(hr0 > hr) break; if(hr0 == hr) if(mn0 >= mn) break; } } if(feof(fi)) {*mysat=0; *a1=*a2=*b1=*b2=0.0;} fclose(fi); free(buf); if(comment) printf("mysat=%d a1=%f b1=%f a2=%f b2=%f\n",*mysat,*a1,*b1,*a2,*b2); } /* eof ---------------------------------------------------------------------------------------------------------------*/