; docformat = 'rst' ; ; NAME: ; shapeMapElevation ; PURPOSE: ; Map the elevation of a body (from its shape model vertices) with respect to a surface reference ;+ ; :Description: ; Map the elevation of a body (from its shape model vertices) with respect to a surface reference ; ; :Categories: ; Shape models ; ; :Returns: An array of elevation maps ; ; :Params: ; model: in, required, type=structure/string ; A shape model structure (see readver) or the path to a shape model ; ; :Keywords: ; abc: in, optional, type=fltarr(3) ; The diameters along the PMI (computed by measureVer if not provided) ; grid: in, optional, type=float, default=0.5 ; The step in longitude/latitude for map creation (in degrees); ; ellipsoid: in, optional, type=string ; If set, write the reference ellipsoid in ver format ; longOff: in, optional, type=float, default=0 ; Move the reference meridian by offSet degrees West. ; ; :Uses: ; measureVer, readVer ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in 2007, B. Carry (ESO):: ; 2016 Dec - B.Carry (OCA) - Full rewrite ; 2017 Dec - B.Carry (OCA) - Removed Dynamical Height ; 2018 Jan - B.Carry (OCA) - Added offset keyword ;- function shapeMapElevation, model, abc=abc, grid=grid, ellipsoid=ellipsoid, longOff=longOff COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Input Syntax Verification -----------------------------------------------------------; if not keyword_set(model) then begin message, /ioError, 'syntax: map = shapeMapElevation( model [, abc=])' return, -1 endif ;--I.2-- Grid Sampling for the Map -----------------------------------------------------------; if not keyword_set(grid) then grid=0.5 nLon = 360. / grid nLat = 180. / grid ;--I.3-- Compute TriAxial Ellipsoid Reference Surface ----------------------------------------; info=measureVer(model) if not keyword_set(abc) then abc=[info.a, info.b, info.c] rMean = product(abc)^(1./3.) ;-Volume equivalent ;--I.4-- Optionnal Longitude Shift of Prime Meridien -----------------------------------------; if not keyword_set(longOff) then longOff=0. ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Create the Reference Ellipsoid -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Convert Vertex Coordinates into Spherical ------------------------------------------; ;--II.1.1-- Shift Vertices to Origin=[0,0,0] xyz = model.vertex.xyz xyz[0,*] -= info.origin[0] xyz[1,*] -= info.origin[1] xyz[2,*] -= info.origin[2] ;--II.1.2-- Spherical to Cartesian sph = cv_coord( from_rect=xyz, /to_sphere ) ;--II.2-- Create an Ellispoid Aligned with XYZ axes ------------------------------------------; ellC = model.vertex.xyz*0. ellC[0,*] = abc[0] * cos(sph[1,*]) * cos(sph[0,*]) ellC[1,*] = abc[1] * cos(sph[1,*]) * sin(sph[0,*]) ellC[2,*] = abc[2] * sin(sph[1,*]) ;--II.3-- Rotate the Ellispoid Along PMI -----------------------------------------------------; t3d, /reset ; t3d, rotate=[0., 0., info.pmi[0].lon ] ; ellR = vert_t3d(ellC) ellR = ellC ; cgPlot, xyz[0,*], xyz[1,*], /Iso, color='Dark Gray';, psym='Filled Circle' ; cgPlot, /OverPlot, ellC[0,*], ellC[1,*], color='Light Gray', psym='Open Circle' ; cgPlot, /OverPlot, ellR[0,*], ellR[1,*], color='Gray';, psym='Filled Circle' ; stop ;--II.4-- Convert Ellipsoid Coordinates into Spherical ---------------------------------------; ell = cv_coord( from_rect=ellR, /to_sphere ) ;--II.5-- Convert Radian into Degree for Mapping ---------------------------------------------; sph[ 0:1, * ] /= !DTOR ell[ 0:1, * ] /= !DTOR ;--II.6-- Move reference Meridian from +180 to 0 ---------------------------------------------; ; sph[0,*] += 180. ; ell[0,*] += 180. ;--II.7-- Introduce an Artificial Noise on Long/Lat to Avoid Pathology with Sph_Scat ---------; nbPoint = n_elements(ell[0,*]) ell[0,*] += randomn(1, nbPoint)/5. ell[1,*] += randomn(2, nbPoint)/5. sph[0,*] += randomn(3, nbPoint)/5. sph[1,*] += randomn(4, nbPoint)/5. ;--II.8-- Export Ellipsoid as VER if Requested -----------------------------------------------; if keyword_set(ellipsoid) then begin ;--II.8.1-- Open File openW, idEll, ellipsoid, /Get_Lun printf, idEll, model.vertex.n, model.facet. n, format='(I,I)' ;--II.8.2-- Write Vertices for k=0, model.vertex.n-1 do $ printf, idEll, ellR[0,k], ellR[1,k], ellR[2,k], format='(F,F,F)' ;--II.8.3-- Write Facets for k=0, model.facet.n-1 do $ printf, idEll, model.facet.list[0,k], model.facet.list[1,k], model.facet.list[2,k], format='(I,I,I)' ;--II.8.4-- Close File free_lun, idEll endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Create the Elevation Maps -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Object Radius ---------------------------------------------------------------------; mapSph = sph_scat( sph[0,*]+longOff, $ ;-Longitude (degree [0,360]) sph[1,*], $ ;-Latitude (degree [-90,90]) sph[2,*], $ ;-Radius [km] nlon=nLon, $ nlat=nLat, $ bounds=[0,-90,360,90.] ) new=congrid(mapSph,nLon,nLat) mapSph=new ;--III.2-- Reference Radius ------------------------------------------------------------------; mapEll = sph_scat( ell[0,*]+longOff, $ ;-Longitude (degree [0,360]) ell[1,*], $ ;-Latitude (degree [-90,90]) ell[2,*], $ ;-Radius [km] nlon=nLon, $ nlat=nLat, $ bounds=[0,-90,360,90.] ) new=congrid(mapEll,nLon,nLat) mapEll=new ;--III.3-- Topography ------------------------------------------------------------------------; mapRel = mapSph - mapEll ; print, sqrt(total( mapRel^2 )/n_elements(mapRel)) ; ; atv, maprel ; ; ; ;longOff=35. ; ; mapSph2 = sph_scat( sph[0,*]+longOff, $ ;-Longitude (degree [0,360]) ; sph[1,*], $ ;-Latitude (degree [-90,90]) ; sph[2,*], $ ;-Radius [km] ; nlon=nLon, $ ; nlat=nLat, $ ; bounds=[0,-90,360,90.]) ; new=congrid(mapSph2,nLon,nLat) ; mapSph2=new ; ; mapEll2 = sph_scat( ell[0,*]+longOff, $ ;-Longitude (degree [0,360]) ; ell[1,*], $ ;-Latitude (degree [-90,90]) ; ell[2,*], $ ;-Radius [km] ; nlon=nLon, $ ; nlat=nLat, $ ; bounds=[0,-90,360,90.] ) ; new=congrid(mapEll2,nLon,nLat) ; mapEll2=new ; mapRel2 = mapSph2 - mapEll2 ; ; diff = maprel-maprel2 ; ; stop ;--III.4-- Dynamical Height ------------------------------------------------------------------; ; mapDyn = rMean*rMean * ( 1./mapEll - 1./mapSph ) ;-- Mauvaise equation. REvoir Thomas et al. 1993 ;help, mapSph, mapEll, mapRel, mapDyn ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Equatorial Cylindrical Projection -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; map = [ [[mapSph]], $ [[mapEll]], $ [[mapRel]] ] return, map end ; stop ; ; help, sph, ellc ;stop ; ;-Checking the Coordinates ; dim = size( vertex ) ; if dim(0) ne 2 then stop ; if dim(1) eq 3 then objetVertex=vertex else objetVertex=transpose(vertex, [1,0]) ; dim = size( vertex ) ; ; ;-Conversion to Spherical Coordinates ; objetSpheric = fltarr( 3, dim(2) ) ; objetSpheric(0,*) = sqrt(total( ObjetVertex^2, 1 )) ; objetSpheric(1,*) = asin( objetVertex(2,*)/ objetSpheric(0,*) ) ; objetSpheric(2,*) = atan( objetVertex(1,*), objetVertex(0,*) ) ; ; ;-Associated tri-axial ellipsoid ; synthVertex = fltarr( 3, dim(2) ) ; synthVertex(0,*) = triAxial(0) * cos(ObjetSpheric(1,*)) * cos(ObjetSpheric(2,*)) ; synthVertex(1,*) = triAxial(1) * cos(ObjetSpheric(1,*)) * sin(ObjetSpheric(2,*)) ; synthVertex(2,*) = triAxial(2) * sin(ObjetSpheric(1,*)) ; synthSpheric = fltarr( 3, dim(2) ) ; synthSpheric(0,*) = sqrt(total( synthVertex^2, 1 )) ; synthSpheric(1,*) = asin( synthVertex(2,*)/ synthSpheric(0,*) ) ; synthSpheric(2,*) = atan( synthVertex(1,*), synthVertex(0,*) ) ; ; ;-Dynamical height ; ; ;-Creating the maps ; mapObjet = SPH_SCAT( $ ; objetSpheric(2,*)/!DTOR+180., $ ;-Longitude (degree [0,360]) ; objetSpheric(1,*)/!DTOR, $ ;-Latitude (degree [-90,90]) ; objetSpheric(0,*), $ ;-Radius ;; gs=spacing, $ ;-Grid ; nlon=360, $ ; nlat=180 ) ; ; mapSynth = SPH_SCAT( $ ; objetSpheric(2,*)/!DTOR+180., $ ;-Longitude (degree [0,360]) ; objetSpheric(1,*)/!DTOR, $ ;-Latitude (degree [-90,90]) ; synthSpheric(0,*), $ ;-Radius ;; gs=spacing, $ ;-Grid ; nlon=360, $ ; nlat=180 ) ;-Grid ; ; mapDyna = SPH_SCAT( $ ; objetSpheric(2,*)/!DTOR+180., $ ;-Longitude (degree [0,360]) ; objetSpheric(1,*)/!DTOR, $ ;-Latitude (degree [-90,90]) ; dynamicHeight, $ ;-Radius ;; gs=spacing, $ ;-Grid ; nlon=360, $ ; nlat=180 ) ;-Grid ; ; ;-Putting the map on a grid starting a 0 deg. long (and not -180) ; dim=size( mapObjet ) ; mapObjet = shift(mapObjet, dim(1)/2., 0.) ; mapSynth = shift(mapSynth, dim(1)/2., 0.) ; mapDyna = shift(mapDyna , dim(1)/2., 0.) ; ; ; ;-Cylindrical projection array ; if keyword_set(equatorial) then $ ; map = [ [[mapObjet]], $ ; [[mapSynth]], $ ; [[mapObjet-mapSynth]], $ ; [[mapDyna]] ] ; ; ;-Orthogonal projection ; if keyword_set(polar) then begin ; ; ;-South pole projection ; geoSouth = MAP_PROJ_INIT('orthographic', $ ; center_latitude=-90., $ ; center_longitude=180.) ; ; southObjet = map_proj_image( mapObjet, $ ; map_structure=geoSouth, $ ; xindex=xSouth, yindex=ySouth, $ ; dimensions=[512,512]) ; southSynth = map_proj_image( mapSynth, $ ; xindex=xSouth, yindex=ySouth, $ ; dimensions=[512,512]) ; southDyna = map_proj_image( mapDyna, $ ; xindex=xSouth, yindex=ySouth, $ ; dimensions=[512,512]) ; ; ; ;-North pole projection ; geoNorth = MAP_PROJ_INIT('orthographic', $ ; center_latitude=90., $ ; center_longitude=180.) ; ; northObjet = map_proj_image( mapObjet, $ ; map_structure=geoNorth, $ ; xindex=xNorth, yindex=yNorth, $ ; dimensions=[512,512]) ; northSynth = map_proj_image( mapSynth, $ ; xindex=xNorth, yindex=yNorth, $ ; dimensions=[512,512]) ; northDyna = map_proj_image( mapDyna, $ ; xindex=xNorth, yindex=yNorth, $ ; dimensions=[512,512]) ; ; ;-Orthographic projection array ; map=[ [[northObjet],[southObjet]], $ ; [[northSynth],[southSynth]], $ ; [[northObjet-northSynth], [southObjet-southSynth]], $ ; [[northDyna],[southDyna]] ] ; ; endif ; ; ; ;-Output ; return, map ; ;end