velocity.py 文件源码

def _read_fmm3d(self, inf):
        inf = open(inf, "r")
        self.nvgrids, self.nvtypes = [int(v) for v in inf.readline().split()[:2]]
        self.v_type_grids = {}
        for (typeID, gridID) in [(ivt, ivg) for ivt in range(1, self.nvtypes+1)
                                            for ivg in range(1, self.nvgrids+1)]:
            if typeID not in self.v_type_grids:
                self.v_type_grids[typeID] = {}
            model = {"typeID": typeID, "gridID": gridID}
            nrho, nlambda, nphi = [int(v) for v in inf.readline().split()[:3]]
            drho, dlambda, dphi = [float(v) for v in inf.readline().split()[:3]]
            rho0, lambda0, phi0 = [float(v) for v in inf.readline().split()[:3]]
            model["grid"] = seispy.geogrid.GeoGrid3D(np.degrees(lambda0),
                                              np.degrees(phi0),
                                              seispy.constants.EARTH_RADIUS - (rho0 + (nrho-1)*drho),
                                              nlambda, nphi, nrho,
                                              np.degrees(dlambda),
                                              np.degrees(dphi),
                                              drho)
            #model["coords"] = coords.SphericalCoordinates((nrho, nlambda, nphi))
            #model["coords"][...] = [[[[rho0 + irho * drho, ?/2 - (lambda0 + ilambda * dlambda), phi0 + iphi * dphi]
            #                           for iphi in range(nphi)]
            #                           for ilambda in range(nlambda)]
            #                           for irho in range(nrho)]
            #model["coords"] = np.flip(model["coords"], axis=1)
            model["data"] = np.empty((nrho, nlambda, nphi))
            model["data"][...] = [[[float(inf.readline().split()[0])
                                    for iphi in range(nphi)]
                                    for ilambda in range(nlambda)]
                                    for irho in range(nrho)]
            model["data"] = np.flip(model["data"], axis=1)
            self.v_type_grids[typeID][gridID] = model