Open sampotter opened 2 years ago
@sampotter I think there is an issue in the "stereo vertices" of your mesh. For example, when you compare the stereographic coordinates of ldem_87s_50mpp
(in meters)
with the stereo coordinates (in km) of the deGerlache clip from earlier in the conversation (https://github.com/sampotter/python-flux/issues/55#issuecomment-1023347306):
one sees that they fit with each other only if x=y and y=-x. Any idea where the issue could come from? The stereo transformation itself seems to be ok, so it might be in the meshing/raster import step.
Once the transformation above is applied, then the illumination/temperature maps look fine (here below an overlap of Tmax in K with the expected PSRs, in yellow - Tmax<160K areas are well aligned with the yellow areas), which also seems to exclude a major problem with the cartesian coordinates (else illumination will be off, etc...)
I would say it is quite possible I did something silly when reading the data in from Mike's GeoTIFF file. I am not sure exactly where the problem is, but I do recall having to screw around with the axes in the image plots to make them match.
Yes, there is some issue with the axes. Already the fact that you have to fix the center of the roi for Gerlache at
xc_roi, yc_roi = 0, -45 # km, stereo SP
(while y=0 for this crater in stereo projection) shows that there is an issue. I don't know if we care about it, but I would like to make sure that the cartesian positions are ok (for direct illumination).
I compared with the results of the rioxarray
library on the same DEM, and the z
are different. We could simply replace the getz
function by this one, but your script requests the elevations z(x,y) one by one (instead of once for an array of x and of y), so that calling that "more complex" function is way too expensive (x1000 run-time per each call). I'm trying to use it as comparison to figure out where is the bug (I think in the indexing within getz, but not sure yet).
Can this be related to the (weird) ordering of the np.array indexes (see result of plt.imshow(dem)
)?
IIRC, Mike's GeoTIFF file contains (buried within it somewhere), a packed 1D array of 64-bit floats corresponding to x coordinates, ditto for y coordinates, and then a packed 2D array of z values. (Maybe the x and y floats are also 2D arrays, but pretty sure it was 1D...)
I am not sure how the 2D array of z values is ordered. Should be either row major or column major. If we can figure that out, I should easily be able to make my getz function consistent with it. Hopefully that will fix the bug.
(I just noticed that xc, yc are not only switched w.r.t. the usual stereo coordinates, but rather xc=-y
and yc=x
- just not visible for deGerlache, since y=0)
IIRC, Mike's GeoTIFF file contains (buried within it somewhere), a packed 1D array of 64-bit floats corresponding to x coordinates, ditto for y coordinates, and then a packed 2D array of z values. (Maybe the x and y floats are also 2D arrays, but pretty sure it was 1D...)
By using
import rioxarray as rio
rds = rio.open_rasterio(tif_path)
x_coords = rds.coords['x'].values
y_coords = rds.coords['y'].values
you should be able to retrieve those coordinates.
(also, FYI, dem = rds.data.values[0]
which should simplify the pipeline)
@steo85it
What has been done with this so far? I would like to push this along