by

Where communities thrive


  • Join over 1.5M+ people
  • Join over 100K+ communities
  • Free without limits
  • Create your own community
People
Repo info
Activity
  • Aug 05 11:07
    MarkWieczorek synchronize #239
  • Aug 05 11:07

    MarkWieczorek on develop

    replace X and Y with shift_orig… (compare)

  • Aug 04 09:47
    MarkWieczorek synchronize #239
  • Aug 04 09:47

    MarkWieczorek on develop

    require matplotlib >=3.3 (compare)

  • Aug 04 09:15
    MarkWieczorek synchronize #239
  • Aug 04 09:15

    MarkWieczorek on develop

    Fix matplotlib variable depreca… (compare)

  • Aug 04 08:34
    MarkWieczorek edited #239
  • Aug 04 08:33
    MarkWieczorek synchronize #239
  • Aug 04 08:33

    MarkWieczorek on develop

    Add release notes for v4.7 (compare)

  • Aug 04 08:21
    MarkWieczorek opened #239
  • Aug 03 21:57

    MarkWieczorek on develop

    Add admittance(), correlation()… add admitance, correlation, and… Add plot_admitcorr method to SH… and 5 more (compare)

  • Aug 03 21:57
    MarkWieczorek closed #238
  • Aug 03 21:57

    MarkWieczorek on develop

    Set xarray.gmt accessor for gri… Merge pull request #237 from Ma… (compare)

  • Aug 03 21:57
    MarkWieczorek closed #237
  • Jul 30 13:04
    MarkWieczorek synchronize #238
  • Jul 25 21:06
    MarkWieczorek opened #238
  • Jul 22 10:34
    megies commented #185
  • Jul 21 20:36
    MarkWieczorek opened #237
  • Jul 21 09:20
    MarkWieczorek closed #236
  • Jul 21 09:08
    wudlike commented #236
Marco Pereira
@ny2292000
@MarkWieczorek I have a very peculiar motivation. I need to properly calculate hyperspherical harmonics to high values of k. Hyperspherical Harmonics have three quantum numbers k,l,m. In principle this would be simple since the hyperspherical harmonics is just a gebenbauer polynomial evaluated at cos(ksi) times the Associate Legendre Polynomial evaluates at cos(theta) times cos(m*phi). The issue, I think, is the normalization factors. They involve large number factorials. I used mpmath. I think I did it rigth...yet... the results tells me that I did it wrong.
Marco Pereira
@ny2292000
So, I suspect there is a lot of room for mistakes and I am sure I am making them. So, in that vein, I would welcome if shtools were to implement Hyperspherical Harmonics. I need the normalization factors properly done to be able to calculate spectra. Considering how much difficult I am having to get this done properly, I wouldn't trust myself to add my work to shtools.
matrixbot
@matrixbot
mark-wieczorek Could you give me a reference as to what these functions are and how they are computed?
mark-wieczorek I probably won't have time to do this myself (I've never had to use them before...), but perhaps someone else would be interested in working on this.
timcurry1
@timcurry1

I am using pyshtools to find the Earth's moon gravity vector to maintain an orbiter vehicle. I have a python script that uses MakeGravGridDH:

degree_order_max = 660 # TEC
degree_order = 330 # TEC
rad, theta, phi, total, pot = gravmag.MakeGravGridDH(
    clm, gm, r0, lmax=degree_order_max, a=constant.a_orbit_moon.value,
    #f=constant.f_moon.value,
    lmax_calc=degree_order,
    omega=constant.omega_moon.value, normal_gravity=0, extend=1) 

My orbiter has a semi major axis around the moon at 88000 km. Can I simply use 88000 km in place of r0? I want to get the magnitude of the gravity vector for each lat/long as the vehicle orbits the moon. I am using GRAIL660 data set.

Mark Wieczorek
@MarkWieczorek
@timcurry1 Hi there. You can do this, but there is one thing to keep in mind: The routine MakeGravGrid was written to compute the gravity on the surface of a planet in a coordinate system rotating with the planet. By specifying omega, the routine will include the centrifugal potential. If you are a spacecraft that is not attached to the planet, then you would want to set omega to zero. This way, the computed gravity will only be from the mass within the Moon, and won't have any rotational components.
If you are really doing orbit computations, then you will probably want to add a gravitational contribution arising from Earth. Also, as orbits are usually eccentric, I'm not sure how you would use this routine in that situation. If the altitude of the spacecraft is really ~80000 km above the surface, I also really doubt that you need to compute the gravity up to degree 600 (maybe ~60 is good enough?), but this is something that you will need to test.
Lastly, you should check out the python classes SHGravCoeffs and SHGravGrid. This might simplify things for you, and will allow you to easily plot maps.
And to answer your question, you should set a to the semimajor axis of the spacecraft, and r0 should be the reference radius of the spherical harmonic coefficients.
timcurry1
@timcurry1

Thanks, Mark. I have been comparing our outputs over the last few weeks. My simulation's gravity magnitude is off by the 5th decimal places from your SHTools at a semi major axis of 3999.9 km. Our values are about 1.0054132 ft/sec^2. I expect that my sim has a different frame. You mentioned that MakeGravGrid uses a fixed frame rotating with the moon. Is there a way to get the Cartesian coordinates and the moon's gravity using SHTools? I am not understanding the flatness attribute while using gravmag.MakeGravGridDH. My thoughts are that my frame does not match yours because of the ecliptic plane, elliptic plane, inertial frame or fixed frame. I am not modeling the earth's influence on the moon's gravity.

I also read in your documentation about errors and anomalies. Can I toggle errors while using the gravmag.MakeGravGridDH?

Thanks for your suggestions.

matrixbot
@matrixbot
mark-wieczorek Hi Tim. A couple quick comments:
mark-wieczorek If the result is off in the 5th decimal place, this is possibly a single vs double precision problem. All computations in SHtools are in double.
mark-wieczorek As for your units of ft/sec^2: What's up with that? 😲
mark-wieczorek And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity of centrifugal potential, MakeGravGrid gives you the gravitational potential from the mass of the Moon in the coordinate system used for the gravity model.
mark-wieczorek * And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity or the centrifugal potential, MakeGravGrid gives you the gravitational potential from the mass of the Moon in the coordinate system used for the gravity model.
mark-wieczorek * And for the reference frame: shtools uses the reference frame of the gravity model. As long as you don't include the normal gravity or the centrifugal potential, MakeGravGrid gives you the gravitational acceleration from the mass of the Moon in the coordinate system used for the gravity model.
matrixbot
@matrixbot
mark-wieczorek Also, you should use the errors of the gravity model only when plotting the power spectrum. The errors aren't supposed to be used to make an map of the uncertainty of the gravity field (for that, you would need the full covariance matrix).
matrixbot
@matrixbot
mark-wieczorek Lastly, for your purposes, just set the flattening to zero. By doing this, you will get a grid of the gravitational acceleration on the surface of a sphere. If you add f /= 0, then you would calculate the acceleration on the surface of a flattened ellipsoid.
timcurry1
@timcurry1

Mark: Our Sim also uses double floating points. I agree with your disappointment in using feet. Our newer sims use metric.

I was looking at your ellipsoid model but I realized that it is not affecting gravity because normal gravity is not used. I assume the reference frame is Z-up through the moon's North pole. Is there a way to extract the reference frame's Cartesian coordinates X, Y and Z? I want to understand how the reference frame is transformed back into earth's ECEF and into the Moon's center inertial (or fixed). I currently have a quaternion between the moon's inertial and fixed and I want to validate it, too.

Thanks for helping me with my validation.

matrixbot
@matrixbot
mark-wieczorek You will have to get an ephemeris for the moon for what you are asking. The standard one is run by JPL. Also be aware that there are different ways to define lunar coordinates. The gravity models use a principal axis frame, but image data sets use a mean sub earth frame. You should look at some of the papers by konopliv et al to see how they define the coordinate system.
Mark Wieczorek
@MarkWieczorek
Hi Everyone. I'm working on adding datasets for the gravity field of Earth to pyshtools, but I'm not sure which ones to include. If you use these, could you let me know what you would like to see? I was thinking of about 4 in total, ranging from models based on spacecraft only data, up to something like egm2008. I looked on the ICGEM website, but there are too many models available...
Ilya Oshchepkov
@ioshchepkov
Hello, Mark. I think, the EGM2008 should definitely be included because it is the de facto standard for now and recommended by IERS Conventions 2010. Although, EGM2020 is coming soon. But we need to be careful with the tidal system. Either it must be specified in the documentation, or in the future we need to add a function to switch from one to another. It is not difficult
Probably, another good choice is GOCO06S - it is a modern satellite-only model. The previous one GOCO05s is recommended for experiments for the International Height Reference Frame (IHRF). As well as XGM2016
Mark Wieczorek
@MarkWieczorek
And what about the EIGEN models?
Someone also suggested GGMplus and WGM2012, but I'm not familiar with either of those...
Ilya Oshchepkov
@ioshchepkov
Also a good choice. All of them are alternatives for the EGM2008. But there is no very big difference between them. In my opinion, as an example, one high-degree model is enough for now. Probably, in the future, we need another repo for working with different data sources, including all ICGEM models
Well, or inside SHTOOLS itself
Mark Wieczorek
@MarkWieczorek
So, in your opinion: EGM2008, GOCO06s, and XGM2016? I don't mind adding more, but at the same time, I think that having too many models might confuse people...
For Topo, I've added Earth2012, and am working on Earth2014 now [the later doesn't have spherical harmonic coefficients for shape models..., so I need both]
Ilya Oshchepkov
@ioshchepkov
Yep, this set seems good enough to me. We can update it in any time, if something much better will arrive
matrixbot
@matrixbot
mark-wieczorek I just talked to Frank Lemoine about this and he suggested a few additional models. This is what I'm going to work on:
  • EGM2008
  • GOCO06s
  • XGM2019e_2159
  • GGM05C
  • EIGEN-6C4
  • EIGEN-GRGS.RL04.MEAN-FIELD
    If this is exessive, let me know!
Mark Wieczorek
@MarkWieczorek
Hi @ioshchepkov, I have a quick question about EGM2008: On the ICGEM web site, there is only a single for that is "tide free". But if you download the original archive from NGA, there are two files for "tide free" and "zero tide". Would you be ok with me only including the single file on ICGEM as a pyshtools dataset, or should I include the two models? ICGEM is probably better for me to use because it takes forever to download anything from the NGA website.
Mark Wieczorek
@MarkWieczorek
Actually, don't worry about it: After 30 minutes I redownloaded their EGM2008 archive and it only had tide free in it. The version I got must have come from somewhere else (and I don't want to investigate where....).
Ilya Oshchepkov
@ioshchepkov
HI, @MarkWieczorek I think ICGEM version is fine. Anyway, the conversion is pretty straight forward. Only the c20 will change. So, there is no point to store both versions
Mark Wieczorek
@MarkWieczorek
Hi @ioshchepkov. I was just looking at the read_icgem_gfc() function that you wrote for pyshtools, and after looking at the gfc files, it looks like there are a few header items that are read, but that aren't returned. This really isn't important, but there is the model name, the normalization, tide_system, and product type. I am just wondering if there is an easy way to return these that won't break your code. For example, I could perhaps add an optional argument like return_header=True that returns would append a dictionary of the full header to the output.
Mark Wieczorek
@MarkWieczorek
...and you might want to take a look at the last couple of commits here:
SHTOOLS/SHTOOLS#231
I added support for reading zip, gzip, and urls in the read_icgem_gfc() function. Let me know if you have any comments.
Mark Wieczorek
@MarkWieczorek
For info, I just merged a big PR onto the develop branch for "datasets". With a single line of code, you can download a spherical harmonic dataset from the original source, and return it to the user as a SHCoeffs class instance. This makes extensive use of pooch (Thanks @leouieda for all the help!), and I have included almost all the research grade topography, gravity, and magnetic field models that you would probably ever want.
I'm happy to add more datasets if I've forgotten any, so just let me know. As part of this, I've also updated some of the spherical harmonic io routines to deal with gz files, zip files, and urls. I also added two new shio routines for binary data and [degree, order, value] formats.
And here is the long list of datasets

Mercury

Topography

  • GTMES150: Topography constructed using laser altimeter and occultation data.

Gravity

  • JGMESS160A : Konolpiv et al. (2020)
  • JGMESS160A_ACCEL : Konolpiv et al. (2020)
  • JGMESS160A_TOPOSIG : Konolpiv et al. (2020)
  • GGMES100V08 : Genova et al. (2019)

Venus

Topography

  • VenusTopo719 : Wieczorek (2015)

Gravity

  • MGNP180U : Konopliv et al. (1999)

The Moon

Topography

  • MoonTopo2600p : Wieczorek (2015)

Gravity

  • GRGM900C : Lemoine et al. (2014)
  • GRGM1200B : Goossens et al. (2020)
  • GRGM1200B_RM1_1E0 : Goossens et al. (2020)
  • GL0900D : Konopliv et al. (2014)
  • GL1500E : Konopliv et al. (2014)

Magnetic Field

  • T2015_449 : Wieczorek (2018), using data from Tsunakawa et al. (2015)

Mars

Topography

  • MarsTopo2600 : Wieczorek (2015)

Gravity

  • GMM3 : Genova et al. (2016)
  • GMM3_RM1_1E0 : Goossens et al. (2017)
  • MRO120D : Konopliv et al. (2016)

Magnetic Field

  • Langlais2019 : Langlais et al. (2019)
  • Morschhauser2014 : Morschhauser et al. (2014)

Vesta

Gravity

  • VESTA20H : Konopliv et al. (2014)

Ceres

Gravity

  • CERES18D : Konopliv et al. (2018)

Earth

Topography

  • Earth2012 : Hirt et al. (2012)
  • Earth2014 : Hirt and Rexer (2014)

Gravity

  • EGM2008 : Pavlis et al. (2012); degree 2190
  • EIGEN_6C4 : Förste et al. (2014); degree 2190
  • GGM05C : Ries et al. (2016); degree 360
  • GOCO06S : Kvas et al. (2019); degree 300, satellite only
  • EIGEN_GRGS_RL04_MEAN_FIELD : Lemoine et al. (2019); degree 300, satellite only
  • XGM2019E : Zingerle et al. (2019); degree 2190

Magnetic field

  • IGRF_13 : Thébault et al. (2015); degree 13
  • SWARM_MLI_2D_0501 : Thébault et al. (2013); degree 133
  • NGDC_720_V3 : Maus (2010); degree 740
  • WDMAM2_800 : Lesur et al. (2016); degree 800
GravMagZhao
@GravMagZhao
Hi @MarkWieczorek According to the regulation of the SHTOOLS, the maximum spherical harmonic degree lmax=nlat/2-1, where nlat is the number of points in latitude. However, as far as I know, the maximum spherical harmonic degree of gravity feld lmax should be equal to nlat. When I do gravity inversion in spherical harmonic domain, I first create spherical harmonics of gravity kernel with degrees 180, and my observation data is 180*360. So in the inversion, the maximum spherical harmonic degree is 89 in this test, and degrees larger than 89 are not involved in the inversion. How can I solve this problem?
matrixbot
@matrixbot
mark-wieczorek Hi GravMagZhao (Gitter) . Have you looked at these instructions that describe the different types of grids that shtools uses?
GravMagZhao
@GravMagZhao
Thanks for your quick reply. Should I use GLQ grid format, rather than DH1 and DH2?
matrixbot
@matrixbot
mark-wieczorek It depends on what you want to do. GLQ is faster by a factor of 2 because there are few latitude bands. But the DH grids are equispaced and easier to work with.
GravMagZhao
@GravMagZhao
My data is equispaced grid, so in this case the maximum spherical harmonic degree is set to be N/2-1.
matrixbot
@matrixbot
mark-wieczorek That is correct. You should read the paper by Driscoll and Healy if you want more details as to why this is the case.
GravMagZhao
@GravMagZhao
Thanks for your suggestion. that is to say, my problem can not be solved by using the SHTOOLS?
matrixbot
@matrixbot
mark-wieczorek If you are using an equispaced grid with nlat = 180, the maximum spherical harmonic degree that can be resolved is 89. That is just how it is.
Mark Wieczorek
@MarkWieczorek
We will be putting out a huge new release in about 2 weeks time. If you want to test this version, please checkout out this PR: SHTOOLS/SHTOOLS#239
Everything is done, and I am just waiting for GMT and pygmt to fix a couple of bugs on their end, and as soon as that is done, I'll merge this into master.
There are lots of new features, of which the addition of datasets is probably the most significant!