Changeset 26367


Ignore:
Timestamp:
07/29/21 16:25:35 (4 years ago)
Author:
jdquinn
Message:

BUG: Corrected CoordTransform; added install script for PROJ 8

Location:
issm/trunk-jpl
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/etc/environment.sh

    r26343 r26367  
    507507if [ -d "${PROJ_ROOT}" ]; then
    508508        export PROJ_ROOT # Used in installation of GDAL
     509        path_append "${PROJ_ROOT}/bin"
    509510        dyld_library_path_append "${PROJ_ROOT}/lib"
    510511        ld_library_path_append "${PROJ_ROOT}/lib"

  • issm/trunk-jpl/src/c/modules/CoordinateSystemTransformx/CoordinateSystemTransformx.cpp

    r25053 r26367  
    44
    55/*Header files*/
     6#include <proj.h>
     7
    68#include "./CoordinateSystemTransformx.h"
    79#include "../../shared/shared.h"
    810#include "../../toolkits/toolkits.h"
    9 #define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
    10 #include <proj_api.h>
    1111
    1212void CoordinateSystemTransformx(double** px_dest,double** py_dest,double* x_src,double* y_src,int size,const char* str_src,const char* str_dst){
    1313
    14 #if !defined(_HAVE_PROJ_)
    15         _error_("proj not installed");
     14/*
     15        NOTE:
     16        - Compilation of this module is fenced in Makefile.am, so we do not
     17        have to check again if _HAVE_PROJ_ is defined
     18        - As of at least PROJ 8, version macros are defined,
     19
     20                https://proj.org/development/reference/macros.html
     21*/
     22#if !defined(PROJ_VERSION_MAJOR) || !PROJ_AT_LEAST_VERSION(8, 0, 0)
     23        _error_("PROJ major version >= 8 required");
    1624#else
     25        /*
     26        Converts an array of longitude (x_src) and array of latitude (y_src) in the
     27        projection described by str_src to the projection described by str_dest.
     28
     29        NOTE:
     30        - API types and calls have been migrated from PROJ 4 to PROJ 6. See SVN
     31        revision history for changes.
     32
     33        Sources:
     34        - https://proj.org/development/migration.html
     35        - https://www.gaia-gis.it/fossil/libspatialite/wiki?name=PROJ.6
     36
     37        TODO:
     38        - Because transformations are done in place, we could save memory and time
     39        by simply passing in x_src and y_src and removing the px_dest and py_dest
     40        parameters entirely?
     41        */
    1742
    1843        /*Allocate output and initialize values as src*/
    … …  
    2045        double* x_dest = xNew<double>(size);
    2146        double* y_dest = xNew<double>(size);
     47
    2248        for(int i=0;i<size;i++){
    2349                x_dest[i] = x_src[i];
    … …  
    2551        }
    2652
    27         /*Create proj.4 projection objects for src and dst*/
    28         projPJ pj_src  = pj_init_plus(str_src);
    29         projPJ pj_dst  = pj_init_plus(str_dst);
    30         if(!pj_src) _error_("Failed to initialize PROJ with source projection \""    <<str_src<<"\"\n");
    31         if(!pj_dst) _error_("Failed to initialize PROJ with destination projection\""<<str_dst<<"\"\n");
     53        PJ *P;
     54        size_t sx = sizeof(double);
     55        size_t sy = sizeof(double);
     56        size_t nx = size;
     57        size_t ny = size;
    3258
    33         /*Perform transformation*/
    34         int p = pj_transform(pj_src,pj_dst,size,1,x_dest,y_dest,NULL);
    35         if(p!=0) _error_("Reprojection failed, PROJ error code: "<<p<<"\n");
     59        P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,str_src,str_dst,NULL);
     60
     61        if ( 0 == P ) {
     62                proj_destroy(P);
     63                _error_("failed to initialize CRS transformation object");
     64        }
     65
     66        int p = proj_trans_generic(P, PJ_FWD, x_dest, sx, nx, y_dest, sy, ny, 0, 0, 0, 0, 0, 0);
     67
     68        if ( 0 == p ){
     69                proj_destroy(P);
     70                _error_("no successful transformations");
     71        }
     72
     73        proj_destroy(P);
    3674
    3775        /*Output : */

  • issm/trunk-jpl/src/m/modules/CoordTransform.m

    r25364 r26367  
    1010%
    1111%   Example of Projections:
    12 %     
     12%
    1313%   lat/lon = '+init=epsg:4326'
    1414%   lat/lon = '+proj=longlat +datum=WGS84'
    15 %   
     15%
    1616%   ll2xy default equivalent (uses with Hugues Ellispoid S)
    1717%   Greenland = '+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.448564109 +units=m +no_defs'
    … …  
    1919%   Bamber Greenland = '+proj=stere +lat_0=90 +lat_ts=71 +lon_0=-39 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs'
    2020%
    21 %   for latitude/longitude, x is longitude and y latitude and the angles are in rad
     21%   For latitude/longitude, x is longitude and y latitude and the angles are in
     22%   rad
    2223%
    23 %   This function will only work of PROJ has been installed in externalpackages
     24%   This function will only work if PROJ has been installed and --with-proj-dir
     25%   option has been set to its location in ISSM configuration
    2426
    2527% Check usage

  • issm/trunk-jpl/src/wrappers/CoordTransform/CoordTransform.cpp

    r25054 r26367  
    33 */
    44
     5/*Header files*/
     6#include <proj.h>
     7
    58#include "./CoordTransform.h"
    6 #define ACCEPT_USE_OF_DEPRECATED_PROJ_API_H
    7 #include <proj_api.h>
    89
    9 void CoordTransformUsage(void){/*{{{*/
    10         _printf_(" type help CoordTransform\n");
    11 }/*}}}*/
    12 WRAPPER(CoordTransform_python){
    13        
    14         /*intermediary: */
    15         double *xin     = NULL;
    16    double *yin     = NULL;
    17    char   *projin  = NULL;
    18    char   *projout = NULL;
    19    int     M,N;
    20    int     test1,test2;
     10/*
     11        NOTE:
     12        - Compilation of this module is fenced in Makefile.am, so we do not
     13        have to check again if _HAVE_PROJ_ is defined
     14        - As of at least PROJ 8, version macros are defined,
    2115
    22         printf("%i %i %i %i\n",NLHS,nlhs,NRHS,nrhs);
     16                https://proj.org/development/reference/macros.html
     17*/
     18#if !defined(PROJ_VERSION_MAJOR)
     19        _error_("PROJ major version >= 8 required");
     20#else
     21        /*
     22        Converts an array of longitude (x_src) and array of latitude (y_src) in the
     23        projection described by str_src to the projection described by str_dest.
    2324
    24         /*Boot module: */
    25         MODULEBOOT();
     25        NOTE:
     26        - API types and calls have been migrated from PROJ 4 to PROJ 6. See SVN
     27        revision history for changes.
    2628
    27         /*checks on arguments: */
    28         CHECKARGUMENTS(NLHS,NRHS,&CoordTransformUsage);
     29        Sources:
     30        - https://proj.org/development/migration.html
     31        - https://www.gaia-gis.it/fossil/libspatialite/wiki?name=PROJ.6
     32        */
    2933
    30         /*Fetch data needed for meshing: */
    31         FetchData(&xin,&M,&N,XIN);
    32    FetchData(&yin,&test1,&test2,YIN);
    33    if(!M*N)     _error_("no coordinate provided");
    34    if(test1!=M) _error_("x and y do not have the same size");
    35    if(test2!=N) _error_("x and y do not have the same size");
    36    FetchData(&projin,PROJIN);
    37    FetchData(&projout,PROJOUT);
     34        void CoordTransformUsage(void){/*{{{*/
     35                _printf_(" type help CoordTransform\n");
     36        }/*}}}*/
     37        WRAPPER(CoordTransform_python){
     38                /*intermediary: */
     39                double *xin     = NULL;
     40                double *yin     = NULL;
     41                char   *projin  = NULL;
     42                char   *projout = NULL;
     43                int     M,N;
     44                int     test1,test2;
    3845
    39         /*Initialize output*/
    40         double* xout = xNew<double>(M*N);
    41         double* yout = xNew<double>(M*N);
     46                /*Boot module: */
     47                MODULEBOOT();
    4248
    43         /*Initialize projections*/
    44         projPJ pj_src = pj_init_plus(projin);
    45         projPJ pj_dst = pj_init_plus(projout);
    46         if(!pj_src) _error_("Failed to initialize PROJ with source projection\n");
    47         if(!pj_dst) _error_("Failed to initialize PROJ with destination projection\n");
     49                /*checks on arguments: */
     50                CHECKARGUMENTS(NLHS,NRHS,&CoordTransformUsage);
    4851
    49         /*Transform coordinates*/
    50         int p = pj_transform(pj_src,pj_dst,M*N,1,xout,yout,NULL);
    51         if(p!=0){
    52                 _error_("PROJ failed with error code: "<<p);
     52                /*Fetch data needed for meshing: */
     53                FetchData(&xin,&M,&N,XIN);
     54                FetchData(&yin,&test1,&test2,YIN);
     55                if(!M*N)     _error_("no coordinate provided");
     56                if(test1!=M) _error_("x and y do not have the same size");
     57                if(test2!=N) _error_("x and y do not have the same size");
     58                FetchData(&projin,PROJIN);
     59                FetchData(&projout,PROJOUT);
     60
     61                /*Calculate size once*/
     62                int size = M*N;
     63
     64                /*Initialize output*/
     65                double* xout = xNew<double>(size);
     66                double* yout = xNew<double>(size);
     67
     68                for(int i=0;i<size;i++){
     69                        xout[i] = xin[i];
     70                        yout[i] = yin[i];
     71                }
     72
     73                PJ *P;
     74                size_t sx = sizeof(double);
     75                size_t sy = sizeof(double);
     76                size_t nx = size;
     77                size_t ny = size;
     78
     79                P = proj_create_crs_to_crs(PJ_DEFAULT_CTX,projin,projout,NULL);
     80
     81                if ( 0 == P ) {
     82                        proj_destroy(P);
     83                        _error_("failed to initialize CRS transformation object");
     84                }
     85
     86                int p = proj_trans_generic(P, PJ_FWD, xout, sx, nx, yout, sy, ny, 0, 0, 0, 0, 0, 0);
     87
     88                if ( 0 == p ){
     89                        proj_destroy(P);
     90                        _error_("no successful transformations");
     91                }
     92
     93                proj_destroy(P);
     94
     95                /*write outputs: */
     96                WriteData(XOUT,xout,M,N);
     97                WriteData(YOUT,yout,M,N);
     98
     99                /*Clean-up and return*/
     100                xDelete<double>(xin);
     101                xDelete<double>(yin);
     102                xDelete<double>(xout);
     103                xDelete<double>(yout);
     104                xDelete<char>(projin);
     105                xDelete<char>(projout);
     106
     107                /*end module: */
     108                MODULEEND();
    53109        }
    54110
    55         /*write outputs: */
    56         WriteData(XOUT,xout,M,N);
    57         WriteData(YOUT,yout,M,N);
    58 
    59         /*Clean-up and return*/
    60    mxFree(projin);
    61    mxFree(projout);
    62         xDelete<double>(xin);
    63         xDelete<double>(yin);
    64         xDelete<double>(xout);
    65         xDelete<double>(yout);
    66         xDelete<char>(projin);
    67         xDelete<char>(projout);
    68 
    69         /*end module: */
    70         MODULEEND();
    71 }
     111#endif

  • issm/trunk-jpl/src/wrappers/CoordTransform/CoordTransform.h

    r25052 r26367  
    1717#endif
    1818
    19 #ifdef _HAVE_JAVASCRIPT_MODULES_
    20 #undef _DO_NOT_LOAD_GLOBALS_ /*only module where this needs to be undefined, so as to
    21                                                            not include IssmComm several times in the javascript Modle construct.*/
    22 #endif
    23 
    2419/*Header files: */
    2520#include "../bindings.h"
     21#include "../../c/main/globals.h"
     22#include "../../c/modules/modules.h"
    2623#include "../../c/shared/shared.h"
    27 #include "../../c/shared/io/io.h"
    2824
    2925#undef __FUNCT__
Note: See TracChangeset for help on using the changeset viewer.