/*---------------------------------------------------------------------------*/ /* 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 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 (3 km, 3712 x 3712 pixels)) image coordinates into low resolution (10 km [actually 9, but called 10 !], 1238 x 1238 pixels) and vice-versa Elevation and azimut of the satellite as seen from the current pixel Time 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. nblig = number of rows (lines) of the image lig = row of the current pixel nbcol = number of columns col = column of the current pixel lat = latitude (decimal degrees) lon = longitude (decimal degrees) I3 denotes the image with the original resolution of 3 km (e.g., I3_Nblig) I3_Nblig = I3_Nbcol = 3712 I10 = same but with the reduced format of 10 km. This format is obtained by a high-pass filtering using a 17x17 kernel, followed by a subsampling (top left of the 3x3 pixels) I10_Nblig = I10_Nbcol = 1238 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 ------------------------------------------------------------------------------*/ /* Some quantities are defined in an external library. They are reported here: */ #define PI 3.141592654 #define R 6371.0 /* Earth radius in km */ #define Io 1367.0 /* solar constant in W/m2 */ #define deg2rad 0.0174533 /* decimal degrees to radian 90°=1.57rd */ #define rad2deg 57.29578 /* converts radians into degrees */ #define Rp 6356.580 /* earth polar radius cf.LRIT/HRIT Global specif. */ #define Re 6378.170 /* earth equatorial radius .......................*/ #define Hs 35785.80 /* altitude of satellite */ #define Lon0 0.0 /* longitude of satellite in degrees */ #define Scanstep 0.00025153/* instrument scan step angle (rd) for 3 lines */ #define I3_NBLIG 3712 /* nblig in 3 km image */ #define I3_NBCOL 3712 /* nbcol ............................. */ #define I10_NBLIG 1238 /* nblig in reduced resolution image (10 km) */ #define I10_NBCOL 1238 /* nbcol ...........................*/ #define Io 1367.0 /* solar constant in W/m2 */ int latlon_to_ligcol(); int ligcol_to_latlon(); int I3_to_latlon(); int latlon_to_I3(); int I10_to_latlon(); int latlon_to_I10(); int I10_to_I3(); int I3_to_I10(); int satelev(); int satazim(); int delta_t(); /*------------------------------------------------------------------------------ Provides the geographical coordinates (degrees, positive NE and negative SW) for any pixel (lig, col)[that is row and column] 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,r2,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 ; 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 ; *col = (int) floor(xr+(nbcol/2)+0.5); *lig = (int) floor(yr+(nblig/2)+0.5); *lig = nblig - *lig +1 ; /* reverse image N-S */ *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; /* reverse image from NW-SE to SE-NW acquisition*/ Y = nblig - lig + 1; X = X - (nbcol/2) - 0.5;/* translate 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 of 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 ; { int ier=0; const double deltaX = Scanstep/3.0*rad2deg; /* 1 pixel scanned /col */ const double deltaY = Scanstep/3.0*rad2deg; /* ............... /line */ ier=ligcol_to_latlon(deltaX,deltaY,I3_NBLIG,I3_NBCOL,lig,col,lat,lon,usable); return ier; } /*----------------------------------------------------------------------------*/ int latlon_to_I3(lat, lon, lig, col, usable) double lat,lon ; int *lig,*col,*usable ; { int ier=0; const double deltaX = Scanstep/3.0*rad2deg; /* 1 pixel scanned /col */ const double deltaY = Scanstep/3.0*rad2deg; /* ............... /line */ ier=latlon_to_ligcol(deltaX,deltaY,I3_NBLIG,I3_NBCOL,lat,lon,lig,col,usable); return ier; } /*----------------------------------------------------------------------------*/ 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; } /*------------------------------------------------------------------------------ Time lag between the beginning of the scanning of the high-resolution image and beginning of scanning the line of interest Lines are scanned 3 at a time. Thus, the time step is 12.5'/3750/3 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 = lig balayées 3x3 Returns 0 if OK, 1 otherwise. */ int delta_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 const double 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; } /* eof -----------------------------------------------------------------------*/