GenericMappingTools / gmt

The Generic Mapping Tools
https://www.generic-mapping-tools.org
Other
862 stars 359 forks source link

3D slice at incorrect position #708

Open zguoch opened 5 years ago

zguoch commented 5 years ago

Description of the problem

image

function preset()
{
    figname=Test_3d
    angle_view=-45/25
    width_fig_x=10
    width_fig_y=10
    width_fig_z=7
    # alpha
    alpha_profile=50
    # range
    xmin=0
    xmax=100
    ymin=0
    ymax=50
    zmin=-20
    zmax=0
    xc=`echo $xmin $xmax | awk '{print $1+($2-$1)/2}'`
    yc=`echo $ymin $ymax | awk '{print $1+($2-$1)/2}'`
    zc=`echo $zmin $zmax | awk '{print $1+($2-$1)/2}'`
    len_z=`echo $zmin $zmax | awk '{print ($2-$1)}'`
}
preset

# plot 
gmt begin $figname pdf,png
    gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Ba -Bza+l"Z(m)" -BwsenZ+gred  -pz$angle_view/$zmin
    echo "$xc $yc zmin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
    gmt basemap -JX$width_fig_x/$width_fig_y -JZ$width_fig_z -R$xmin/$xmax/$ymin/$ymax/$zmin/$zmax -Bwsenz+ggray  -pz$angle_view/$zmax
    echo "$xc $yc zmax plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a190 -Dj0c/0c
    gmt basemap -JX$width_fig_y/$width_fig_z -JZ$width_fig_x -R$ymin/$ymax/$zmin/$zmax/$xmin/$xmax -Ba -BwSenz+glightblue@$alpha_profile -Bx+l"Y(m)" --MAP_FRAME_PEN=1,black@0 --MAP_ANNOT_OFFSET=-0.2   -px$angle_view/$xminn
    echo "$yc $zc ymin plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a45 -Dj0c/0c
    gmt basemap -JX$width_fig_x/$width_fig_z -JZ$width_fig_y -R$xmin/$xmax/$zmin/$zmax/$ymin/$ymax -Ba -BwSenz+glightgreen@$alpha_profile -Bx+l"X(m)" --MAP_FRAME_PEN=1,black@0 --MAP_ANNOT_OFFSET=-0.2  -py$angle_view/$ymax
    echo "$xc $zc xminn plane" | gmt pstext -JZ -p -F+f20p,Helvetica-Bold,blue=thinner,white+jCM+a-45 -Dj0c/0c
gmt end
open $figname.pdf
rm tmp* gmt.conf gmt.history 

System information

welcome[bot] commented 5 years ago

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. We appreciate that you took the time to contribute!

Please make sure you read our Contributing Guide and abide by our Code of Conduct.

zguoch commented 5 years ago

I will take time to read the Contributing Guid and Code of Conduct, and will make a "pull request" for this issue. Here I attach the key part of code to solve this issue.

GMT_LOCAL int map_init_three_D (struct GMT_CTRL *GMT) {
//...
//...
if (easy) {
        double xx[4], yy[4];

        xx[0] = xx[3] = GMT->current.proj.rect[XLO]; xx[1] = xx[2] = GMT->current.proj.rect[XHI];
        yy[0] = yy[1] = GMT->current.proj.rect[YLO]; yy[2] = yy[3] = GMT->current.proj.rect[YHI];
        // which plane in 3d projection
        int plane = GMT->current.proj.z_project.view_plane;
        // the offset coefficenece of transformation matrix: e and e should be keep consistent with the case of GMT_Z
        switch (plane % 3) {
            case GMT_X:
                // -Rymin/ymax/zmin/zmax/xmin/xmax
                xx[0] = xx[3] = GMT->current.proj.zmin; xx[1] = xx[2] = GMT->current.proj.zmax;
                yy[0] = yy[1] = GMT->current.proj.rect[XLO]; yy[2] = yy[3] = GMT->current.proj.rect[XHI];
                break;
            case GMT_Y:
                // -Rxmin/xmax/zmin/zmax/ymin/ymax
                xx[0] = xx[3] = GMT->current.proj.rect[XLO]; xx[1] = xx[2] = GMT->current.proj.rect[XHI];
                yy[0] = yy[1] = GMT->current.proj.zmin; yy[2] = yy[3] = GMT->current.proj.zmax;
                break;
            case GMT_Z:
                // -Rxmin/xmax/ymin/ymax/zmin/zmax
                xx[0] = xx[3] = GMT->current.proj.rect[XLO]; xx[1] = xx[2] = GMT->current.proj.rect[XHI];
                yy[0] = yy[1] = GMT->current.proj.rect[YLO]; yy[2] = yy[3] = GMT->current.proj.rect[YHI];
                break;
        }
        for (i = 0; i < 4; i++) {
            gmt_xy_to_geo (GMT, &GMT->current.proj.z_project.corner_x[i], &GMT->current.proj.z_project.corner_y[i], xx[i], yy[i]);
            gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZLO]), &x, &y);
            GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x);
            GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x);
            GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y);
            GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y);
            gmt_xyz_to_xy (GMT, xx[i], yy[i], gmt_z_to_zz(GMT, GMT->common.R.wesn[ZHI]), &x, &y);
            GMT->current.proj.z_project.xmin = MIN (GMT->current.proj.z_project.xmin, x);
            GMT->current.proj.z_project.xmax = MAX (GMT->current.proj.z_project.xmax, x);
            GMT->current.proj.z_project.ymin = MIN (GMT->current.proj.z_project.ymin, y);
            GMT->current.proj.z_project.ymax = MAX (GMT->current.proj.z_project.ymax, y);
        }
    }
//...
}
PaulWessel commented 5 years ago

What you say, @remkos ?

remkos commented 5 years ago

This looks like a duplicate of issue #529. I'll put it on me to implement and test the code by @CosmicScholar Thank you for your contribution!

joa-quim commented 1 year ago

It's a complicated test (for me) to reproduce. Is it still buggy (there changes in this regard) or can we close it?

PaulWessel commented 1 year ago

This problem seems to have taken a turn for the worse. Running it today I get

Test_3d