Skip to content
GitLab
Commit 3cc3199e authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added function cdo_lcc_to_lonlat().

parent 06e9966a
Hide whitespace changes
Inline Side-by-side
src/Makefile.am
View file @ 3cc3199e
......@@ -66,6 +66,7 @@ libcdo_la_SOURCES = \
gradsdeslib.h \
grid.c \
grid.h \
grid_proj.c \
grid_area.c \
grid_define.c \
grid_gme.c \
......
src/Makefile.in
View file @ 3cc3199e
......@@ -118,16 +118,17 @@ am_libcdo_la_OBJECTS = libcdo_la-array.lo libcdo_la-cdo_pthread.lo \
libcdo_la-field.lo libcdo_la-field2.lo libcdo_la-fieldc.lo \
libcdo_la-fieldmem.lo libcdo_la-fieldmer.lo \
libcdo_la-fieldzon.lo libcdo_la-gradsdeslib.lo \
libcdo_la-grid.lo libcdo_la-grid_area.lo \
libcdo_la-grid_define.lo libcdo_la-grid_gme.lo \
libcdo_la-grid_rot.lo libcdo_la-grid_from_name.lo \
libcdo_la-grid_read.lo libcdo_la-grid_read_pingo.lo \
libcdo_la-grid_print.lo libcdo_la-gridreference.lo \
libcdo_la-griddes.lo libcdo_la-griddes_h5.lo \
libcdo_la-griddes_nc.lo libcdo_la-hetaeta.lo \
libcdo_la-institution.lo libcdo_la-interpol.lo \
libcdo_la-job.lo libcdo_la-juldate.lo libcdo_la-grid_search.lo \
libcdo_la-listarray.lo libcdo_la-list.lo libcdo_la-listbuf.lo \
libcdo_la-grid.lo libcdo_la-grid_proj.lo \
libcdo_la-grid_area.lo libcdo_la-grid_define.lo \
libcdo_la-grid_gme.lo libcdo_la-grid_rot.lo \
libcdo_la-grid_from_name.lo libcdo_la-grid_read.lo \
libcdo_la-grid_read_pingo.lo libcdo_la-grid_print.lo \
libcdo_la-gridreference.lo libcdo_la-griddes.lo \
libcdo_la-griddes_h5.lo libcdo_la-griddes_nc.lo \
libcdo_la-hetaeta.lo libcdo_la-institution.lo \
libcdo_la-interpol.lo libcdo_la-job.lo libcdo_la-juldate.lo \
libcdo_la-grid_search.lo libcdo_la-listarray.lo \
libcdo_la-list.lo libcdo_la-listbuf.lo \
libcdo_la-merge_sort2.lo libcdo_la-modules.lo \
libcdo_la-namelist.lo libcdo_la-namelist_parser.lo \
libcdo_la-normal.lo libcdo_la-nth_element.lo \
......@@ -522,15 +523,16 @@ libcdo_la_SOURCES = array.h array.c cdo_int.h compare.h cdo_pthread.c \
expr_fun.c expr_fun.h expr_lex.c expr_yacc.c expr_yacc.h \
features.c field.c field.h field2.c fieldc.c fieldmem.c \
fieldmer.c fieldzon.c functs.h gradsdeslib.c gradsdeslib.h \
grid.c grid.h grid_area.c grid_define.c grid_gme.c grid_rot.c \
grid_from_name.c grid_read.c grid_read_pingo.c grid_print.c \
gridreference.c griddes.c griddes.h griddes_h5.c griddes_nc.c \
hetaeta.c hetaeta.h institution.c interpol.c interpol.h job.c \
juldate.c grid_search.c grid_search.h listarray.c listarray.h \
list.c list.h listbuf.c listbuf.h merge_sort2.c merge_sort2.h \
modules.c modules.h namelist.c namelist.h namelist_parser.c \
normal.c nth_element.c nth_element.h operator_help.h par_io.c \
par_io.h parse_literal.c percentiles_hist.c percentiles_hist.h \
grid.c grid.h grid_proj.c grid_area.c grid_define.c grid_gme.c \
grid_rot.c grid_from_name.c grid_read.c grid_read_pingo.c \
grid_print.c gridreference.c griddes.c griddes.h griddes_h5.c \
griddes_nc.c hetaeta.c hetaeta.h institution.c interpol.c \
interpol.h job.c juldate.c grid_search.c grid_search.h \
listarray.c listarray.h list.c list.h listbuf.c listbuf.h \
merge_sort2.c merge_sort2.h modules.c modules.h namelist.c \
namelist.h namelist_parser.c normal.c nth_element.c \
nth_element.h operator_help.h par_io.c par_io.h \
parse_literal.c percentiles_hist.c percentiles_hist.h \
percentiles.c percentiles.h pipe.c pipe.h pmlist.c pmlist.h \
sellist.c sellist.h pragma_omp_atomic_update.h printinfo.h \
process.c process.h pstream.c pstream.h pstream_write.h \
......@@ -1051,6 +1053,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_from_name.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_gme.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_print.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_proj.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_read.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_read_pingo.Plo@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/libcdo_la-grid_rot.Plo@am__quote@
......@@ -1390,6 +1393,13 @@ libcdo_la-grid.lo: grid.c
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libcdo_la-grid.lo `test -f 'grid.c' || echo '$(srcdir)/'`grid.c
 
libcdo_la-grid_proj.lo: grid_proj.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libcdo_la-grid_proj.lo -MD -MP -MF $(DEPDIR)/libcdo_la-grid_proj.Tpo -c -o libcdo_la-grid_proj.lo `test -f 'grid_proj.c' || echo '$(srcdir)/'`grid_proj.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libcdo_la-grid_proj.Tpo $(DEPDIR)/libcdo_la-grid_proj.Plo
@AMDEP_TRUE@@am__fastdepCC_FALSE@ $(AM_V_CC)source='grid_proj.c' object='libcdo_la-grid_proj.lo' libtool=yes @AMDEPBACKSLASH@
@AMDEP_TRUE@@am__fastdepCC_FALSE@ DEPDIR=$(DEPDIR) $(CCDEPMODE) $(depcomp) @AMDEPBACKSLASH@
@am__fastdepCC_FALSE@ $(AM_V_CC@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -c -o libcdo_la-grid_proj.lo `test -f 'grid_proj.c' || echo '$(srcdir)/'`grid_proj.c
libcdo_la-grid_area.lo: grid_area.c
@am__fastdepCC_TRUE@ $(AM_V_CC)$(LIBTOOL) $(AM_V_lt) --tag=CC $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CC) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libcdo_la_CPPFLAGS) $(CPPFLAGS) $(AM_CFLAGS) $(CFLAGS) -MT libcdo_la-grid_area.lo -MD -MP -MF $(DEPDIR)/libcdo_la-grid_area.Tpo -c -o libcdo_la-grid_area.lo `test -f 'grid_area.c' || echo '$(srcdir)/'`grid_area.c
@am__fastdepCC_TRUE@ $(AM_V_at)$(am__mv) $(DEPDIR)/libcdo_la-grid_area.Tpo $(DEPDIR)/libcdo_la-grid_area.Plo
......
src/grid.c
View file @ 3cc3199e
......@@ -675,151 +675,6 @@ void grid_inq_param_laea(int gridID, double *a, double *lon_0, double *lat_0, do
}
}
#if defined(HAVE_LIBPROJ)
static
bool grid_inq_param_lcc(int gridID, double *a, double *rf, double *lon_0, double *lat_0, double *lat_1, double *lat_2, double *x_0, double *y_0)
{
bool status = false;
double xlon_0 = PARAM_MISSVAL, ylat_0 = PARAM_MISSVAL;
*a = PARAM_MISSVAL; *rf = PARAM_MISSVAL; *lon_0 = PARAM_MISSVAL; *lat_0 = PARAM_MISSVAL;
*lat_1 = PARAM_MISSVAL, *lat_2 = PARAM_MISSVAL;
*x_0 = PARAM_MISSVAL, *y_0 = PARAM_MISSVAL;
int gridtype = gridInqType(gridID);
if ( gridtype == GRID_PROJECTION )
{
const char *projection = "lambert_conformal_conic";
char mapping[CDI_MAX_NAME]; mapping[0] = 0;
cdiGridInqKeyStr(gridID, CDI_KEY_MAPNAME, CDI_MAX_NAME, mapping);
if ( mapping[0] && strcmp(mapping, projection) == 0 )
{
int atttype, attlen;
char attname[CDI_MAX_NAME+1];
int natts;
cdiInqNatts(gridID, CDI_GLOBAL, &natts);
for ( int iatt = 0; iatt < natts; ++iatt )
{
cdiInqAtt(gridID, CDI_GLOBAL, iatt, attname, &atttype, &attlen);
if ( attlen > 2 ) continue;
double attflt[2];
if ( cdiInqAttConvertedToFloat(gridID, atttype, attname, attlen, attflt) )
{
if ( strcmp(attname, "earth_radius") == 0 ) *a = attflt[0];
else if ( strcmp(attname, "inverse_flattening") == 0 ) *rf = attflt[0];
else if ( strcmp(attname, "longitude_of_central_meridian") == 0 ) *lon_0 = attflt[0];
else if ( strcmp(attname, "latitude_of_projection_origin") == 0 ) *lat_0 = attflt[0];
else if ( strcmp(attname, "false_easting") == 0 ) *x_0 = attflt[0];
else if ( strcmp(attname, "false_northing") == 0 ) *y_0 = attflt[0];
else if ( strcmp(attname, "longitudeOfFirstGridPointInDegrees") == 0 ) xlon_0 = attflt[0];
else if ( strcmp(attname, "latitudeOfFirstGridPointInDegrees") == 0 ) ylat_0 = attflt[0];
else if ( strcmp(attname, "standard_parallel") == 0 )
{
*lat_1 = attflt[0];
*lat_2 = (attlen == 2) ? attflt[1] : attflt[0];
}
}
}
status = true;
if ( IS_EQUAL(*lon_0,PARAM_MISSVAL) ) { status = false; cdoWarning("%s mapping parameter %s missing!", projection, "longitude_of_central_meridian"); }
if ( IS_EQUAL(*lat_0,PARAM_MISSVAL) ) { status = false; cdoWarning("%s mapping parameter %s missing!", projection, "latitude_of_central_meridian"); }
if ( IS_EQUAL(*lat_1,PARAM_MISSVAL) ) { status = false; cdoWarning("%s mapping parameter %s missing!", projection, "standard_parallel"); }
if ( status && IS_EQUAL(*x_0,PARAM_MISSVAL) && IS_EQUAL(*y_0,PARAM_MISSVAL) && IS_NOT_EQUAL(xlon_0,PARAM_MISSVAL) && IS_NOT_EQUAL(ylat_0,PARAM_MISSVAL) )
{ status = false; cdoWarning("%s mapping parameter %s missing!", projection, "false_easting and false_northing"); }
}
else
cdoWarning("%s mapping parameter missing!", projection);
}
return status;
}
#endif
int proj_lonlat_to_lcc(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, int nvals, double *xvals, double *yvals)
{
int status = 0;
#if defined(HAVE_LIBPROJ)
char *params[20];
int nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( IS_NOT_EQUAL(a, missval) && a > 0 ) params[nbpar++] = gen_param("a=%g", a);
if ( IS_NOT_EQUAL(rf, missval) && rf > 0 ) params[nbpar++] = gen_param("rf=%g", rf);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
params[nbpar++] = gen_param("units=m");
// params[nbpar++] = gen_param("no_defs");
if ( cdoVerbose )
for ( int i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
projPJ proj = pj_init(nbpar, &params[0]);
if ( !proj ) status = -1;
for ( int i = 0; i < nbpar; ++i ) Free(params[i]);
/* proj->over = 1; */ /* allow longitude > 180 */
if ( status == 0 )
{
projUV p;
for ( int i = 0; i < nvals; i++ )
{
p.u = xvals[i]*DEG_TO_RAD;
p.v = yvals[i]*DEG_TO_RAD;
p = pj_fwd(p, proj);
xvals[i] = p.u;
yvals[i] = p.v;
}
pj_free(proj);
}
#else
status = -1;
#endif
if ( status == -1 )
for ( int i = 0; i < nvals; i++ )
{
xvals[i] = missval;
yvals[i] = missval;
}
return status;
}
static
void lonlat_to_lcc(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, int nvals, double *xvals, double *yvals)
{
int status = proj_lonlat_to_lcc(missval, lon_0, lat_0, lat_1, lat_2, a, rf, nvals, xvals, yvals);
#if defined(HAVE_LIBPROJ)
if ( status == -1 ) cdoAbort("proj error: %s", pj_strerrno(pj_errno));
#else
if ( status == -1 ) cdoAbort("proj4 support not compiled in!");
#endif
}
int cdo_lonlat_to_lcc(int gridID, int nvals, double *xvals, double *yvals)
{
double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;
gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
lonlat_to_lcc(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, nvals, xvals, yvals);
return 0;
}
static
void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
{
......@@ -866,57 +721,6 @@ void laea_to_geo(int gridID, int gridsize, double *xvals, double *yvals)
#endif
}
int lcc_to_lonlat(int gridID, int gridsize, double *xvals, double *yvals)
{
bool errstatus = false;
#if defined(HAVE_LIBPROJ)
char *params[20];
double a, rf, lon_0, lat_0, lat_1, lat_2, x_0, y_0;
bool status = grid_inq_param_lcc(gridID, &a, &rf, &lon_0, &lat_0, &lat_1, &lat_2, &x_0, &y_0);
if ( status == false ) cdoAbort("mapping parameter missing!");
int nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( IS_NOT_EQUAL(a,PARAM_MISSVAL) && a > 0 ) params[nbpar++] = gen_param("a=%g", a);
if ( IS_NOT_EQUAL(rf,PARAM_MISSVAL) && rf > 0 ) params[nbpar++] = gen_param("rf=%g", rf);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
if ( IS_NOT_EQUAL(x_0,PARAM_MISSVAL) ) params[nbpar++] = gen_param("x_0=%g", x_0);
if ( IS_NOT_EQUAL(y_0,PARAM_MISSVAL) ) params[nbpar++] = gen_param("y_0=%g", y_0);
if ( cdoVerbose )
for ( int i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
projPJ proj = pj_init(nbpar, &params[0]);
if ( !proj ) cdoAbort("proj error: %s", pj_strerrno(pj_errno));
for ( int i = 0; i < nbpar; ++i ) Free(params[i]);
/* proj->over = 1; */ /* allow longitude > 180 */
projUV p;
for ( int i = 0; i < gridsize; i++ )
{
p.u = xvals[i];
p.v = yvals[i];
p = pj_inv(p, proj);
xvals[i] = p.u*RAD_TO_DEG;
yvals[i] = p.v*RAD_TO_DEG;
}
pj_free(proj);
#else
errstatus = true;
cdoAbort("proj4 support not compiled in!");
#endif
return !errstatus;
}
static
void proj_to_geo(char *proj4param, int gridsize, double *xvals, double *yvals)
{
......@@ -1411,7 +1215,7 @@ int gridToCurvilinear(int gridID1, int lbounds)
if ( lproj_sinu ) sinu_to_geo(gridsize, xvals2D, yvals2D);
else if ( lproj_laea ) laea_to_geo(gridID1, gridsize, xvals2D, yvals2D);
else if ( lproj_lcc ) lcc_to_lonlat(gridID1, gridsize, xvals2D, yvals2D);
else if ( lproj_lcc ) cdo_lcc_to_lonlat(gridID1, gridsize, xvals2D, yvals2D);
else if ( lproj4 ) proj_to_geo(proj4param, gridsize, xvals2D, yvals2D);
}
......@@ -1498,7 +1302,7 @@ int gridToCurvilinear(int gridID1, int lbounds)
if ( lproj_sinu ) sinu_to_geo(4*gridsize, xbounds2D, ybounds2D);
else if ( lproj_laea ) laea_to_geo(gridID1, 4*gridsize, xbounds2D, ybounds2D);
else if ( lproj_lcc ) lcc_to_lonlat(gridID1, 4*gridsize, xbounds2D, ybounds2D);
else if ( lproj_lcc ) cdo_lcc_to_lonlat(gridID1, 4*gridsize, xbounds2D, ybounds2D);
else if ( lproj4 ) proj_to_geo(proj4param, 4*gridsize, xbounds2D, ybounds2D);
}
else
......
src/grid.h
View file @ 3cc3199e
......@@ -103,11 +103,16 @@ int cdo_define_sample_grid(int gridID, int sampleFactor);
// Define a sub-grid of another grid
int cdo_define_subgrid_grid(int gridSrcID, int subI0, int subI1, int subJ0, int subJ1);
int cdo_lonlat_to_lcc(int gridID, size_t nvals, double *xvals, double *yvals);
int cdo_lcc_to_lonlat(int gridID, size_t nvals, double *xvals, double *yvals);
#ifdef __cplusplus
extern "C" {
#endif
int proj_lonlat_to_lcc(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, int nvals, double *xvals, double *yvals);
double a, double rf, size_t nvals, double *xvals, double *yvals);
#if defined (__cplusplus)
}
#endif
......
src/grid_proj.c 0 → 100644
View file @ 3cc3199e
/*
This file is part of CDO. CDO is a collection of Operators to
manipulate and analyse Climate model Data.
Copyright (C) 2003-2017 Uwe Schulzweida, <uwe.schulzweida AT mpimet.mpg.de>
See COPYING file for copying and redistribution conditions.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
*/
#if defined(HAVE_CONFIG_H)
#include "config.h"
#endif
#if defined(_OPENMP)
#include <omp.h>
#endif
#if defined(HAVE_LIBPROJ)
#include "proj_api.h"
#endif
#include <stdio.h>
#include <stdarg.h> /* va_list */
#include <cdi.h>
#include "cdo_int.h"
#include "grid.h"
static
char *gen_param(const char *fmt, ...)
{
va_list args;
char str[256];
va_start(args, fmt);
int len = vsprintf(str, fmt, args);
va_end(args);
len++;
char *rstr = (char*) Malloc(len*sizeof(char));
memcpy(rstr, str, len*sizeof(char));
return rstr;
}
int proj_lonlat_to_lcc(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, size_t nvals, double *xvals, double *yvals)
{
int status = 0;
#if defined(HAVE_LIBPROJ)
char *params[20];
int nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( IS_NOT_EQUAL(a, missval) && a > 0 ) params[nbpar++] = gen_param("a=%g", a);
if ( IS_NOT_EQUAL(rf, missval) && rf > 0 ) params[nbpar++] = gen_param("rf=%g", rf);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
params[nbpar++] = gen_param("units=m");
// params[nbpar++] = gen_param("no_defs");
if ( cdoVerbose )
for ( int i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
projPJ proj = pj_init(nbpar, &params[0]);
if ( !proj ) status = -1;
for ( int i = 0; i < nbpar; ++i ) Free(params[i]);
/* proj->over = 1; */ /* allow longitude > 180 */
if ( status == 0 )
{
projUV p;
for ( size_t i = 0; i < nvals; i++ )
{
p.u = xvals[i]*DEG_TO_RAD;
p.v = yvals[i]*DEG_TO_RAD;
p = pj_fwd(p, proj);
xvals[i] = p.u;
yvals[i] = p.v;
}
pj_free(proj);
}
#else
status = -1;
#endif
if ( status == -1 )
for ( size_t i = 0; i < nvals; i++ )
{
xvals[i] = missval;
yvals[i] = missval;
}
return status;
}
static
void lonlat_to_lcc(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, size_t nvals, double *xvals, double *yvals)
{
int status = proj_lonlat_to_lcc(missval, lon_0, lat_0, lat_1, lat_2, a, rf, nvals, xvals, yvals);
#if defined(HAVE_LIBPROJ)
if ( status == -1 ) cdoAbort("proj error: %s", pj_strerrno(pj_errno));
#else
if ( status == -1 ) cdoAbort("proj4 support not compiled in!");
#endif
}
int cdo_lonlat_to_lcc(int gridID, size_t nvals, double *xvals, double *yvals)
{
double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;
gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
lonlat_to_lcc(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, nvals, xvals, yvals);
return 0;
}
int proj_lcc_to_lonlat(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, double x_0, double y_0, size_t nvals, double *xvals, double *yvals)
{
int status = 0;
#if defined(HAVE_LIBPROJ)
char *params[20];
int nbpar = 0;
params[nbpar++] = gen_param("proj=lcc");
if ( IS_NOT_EQUAL(a, grid_missval) && a > 0 ) params[nbpar++] = gen_param("a=%g", a);
if ( IS_NOT_EQUAL(rf, grid_missval) && rf > 0 ) params[nbpar++] = gen_param("rf=%g", rf);
params[nbpar++] = gen_param("lon_0=%g", lon_0);
params[nbpar++] = gen_param("lat_0=%g", lat_0);
params[nbpar++] = gen_param("lat_1=%g", lat_1);
params[nbpar++] = gen_param("lat_2=%g", lat_2);
if ( IS_NOT_EQUAL(x_0, grid_missval) ) params[nbpar++] = gen_param("x_0=%g", x_0);
if ( IS_NOT_EQUAL(y_0, grid_missval) ) params[nbpar++] = gen_param("y_0=%g", y_0);
if ( cdoVerbose )
for ( int i = 0; i < nbpar; ++i )
cdoPrint("Proj.param[%d] = %s", i+1, params[i]);
projPJ proj = pj_init(nbpar, &params[0]);
if ( !proj ) status = -1;
for ( int i = 0; i < nbpar; ++i ) Free(params[i]);
/* proj->over = 1; */ /* allow longitude > 180 */
if ( status == 0 )
{
projUV p;
for ( size_t i = 0; i < nvals; i++ )
{
p.u = xvals[i];
p.v = yvals[i];
p = pj_inv(p, proj);
xvals[i] = p.u*RAD_TO_DEG;
yvals[i] = p.v*RAD_TO_DEG;
}
pj_free(proj);
}
#else
status = -1;
#endif
if ( status == -1 )
for ( size_t i = 0; i < nvals; i++ )
{
xvals[i] = missval;
yvals[i] = missval;
}
return status;
}
static
void lcc_to_lonlat(double missval, double lon_0, double lat_0, double lat_1, double lat_2,
double a, double rf, double x_0, double y_0, size_t nvals, double *xvals, double *yvals)
{
int status = proj_lcc_to_lonlat(missval, lon_0, lat_0, lat_1, lat_2, a, rf, x_0, y_0, nvals, xvals, yvals);
#if defined(HAVE_LIBPROJ)
if ( status == -1 ) cdoAbort("proj error: %s", pj_strerrno(pj_errno));
#else
if ( status == -1 ) cdoAbort("proj4 support not compiled in!");
#endif
}
int cdo_lcc_to_lonlat(int gridID, size_t nvals, double *xvals, double *yvals)
{
const char *projection = "lambert_conformal_conic";
double lon_0, lat_0, lat_1, lat_2, a, rf, xval_0, yval_0, x_0, y_0;
gridInqParamLCC(gridID, grid_missval, &lon_0, &lat_0, &lat_1, &lat_2, &a, &rf, &xval_0, &yval_0, &x_0, &y_0);
int status = 0;
if ( !status && IS_EQUAL(lon_0, grid_missval) ) { status = 1; cdoWarning("%s mapping parameter %s missing!", projection, "longitude_of_central_meridian"); }
if ( !status && IS_EQUAL(lat_0, grid_missval) ) { status = 1; cdoWarning("%s mapping parameter %s missing!", projection, "latitude_of_central_meridian"); }
if ( !status && IS_EQUAL(lat_1, grid_missval) ) { status = 1; cdoWarning("%s mapping parameter %s missing!", projection, "standard_parallel"); }
if ( !status && IS_EQUAL(x_0, grid_missval) && IS_EQUAL(y_0, grid_missval) && IS_NOT_EQUAL(xval_0, grid_missval) && IS_NOT_EQUAL(yval_0, grid_missval) )
{
#if defined(HAVE_LIBPROJ)
x_0 = xval_0; y_0 = yval_0;
lonlat_to_lcc(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, 1, &x_0, &y_0);
x_0 = -x_0; y_0 = -y_0;
#else
status = 1; cdoWarning("%s mapping parameter %s missing!", projection, "false_easting and false_northing");
#endif
}
if ( status ) cdoAbort("%s mapping parameter missing!", projection);
lcc_to_lonlat(grid_missval, lon_0, lat_0, lat_1, lat_2, a, rf, x_0, y_0, nvals, xvals, yvals);
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Imprint | Privacy Policy