From 430c1eeb66d02c14af88083e0c5ad7c4a72e7c20 Mon Sep 17 00:00:00 2001 From: rsiddans Date: Thu, 28 Nov 2024 16:44:22 +0000 Subject: [PATCH] Fixes to netCDF4 and hdf5 (#1926) Co-authored-by: Richard Siddans Co-authored-by: Giloo <33936193+GillesDuvert@users.noreply.github.com> --- src/hdf5_fun.cpp | 30 +++++--- src/libinit_cl.cpp | 4 +- src/ncdf_att_cl.cpp | 126 ++++++++++++++++++++++++++++++-- src/ncdf_cl.cpp | 4 ++ src/ncdf_var_cl.cpp | 130 +++++++++++++++++++++++++++++++--- src/pro/map/map_proj_init.pro | 3 +- 6 files changed, 271 insertions(+), 26 deletions(-) diff --git a/src/hdf5_fun.cpp b/src/hdf5_fun.cpp index 314e837d4..8e7439ba0 100644 --- a/src/hdf5_fun.cpp +++ b/src/hdf5_fun.cpp @@ -155,10 +155,10 @@ namespace lib { // { H5T_NATIVE_LDOUBLE }, //not until LDOUBLE Is handled everywhere!! /* GDL_DOUBLE */ - { H5T_NATIVE_DOUBLE }, + { H5T_NATIVE_DOUBLE,H5T_IEEE_F64BE,H5T_IEEE_F64LE,H5T_INTEL_F64,H5T_MIPS_F64,H5T_ALPHA_F64}, /* GDL_FLOAT */ - { H5T_NATIVE_FLOAT }, + { H5T_NATIVE_FLOAT,H5T_ALPHA_F32,H5T_IEEE_F32BE, H5T_IEEE_F32LE,H5T_INTEL_F32,H5T_MIPS_F32}, /* GDL_ULONG64 */ { H5T_NATIVE_ULLONG, H5T_ALPHA_U64, H5T_INTEL_U64, H5T_MIPS_U64, @@ -166,12 +166,12 @@ namespace lib { H5T_STD_U64BE, H5T_STD_U64LE }, /* GDL_LONG64 */ - { H5T_NATIVE_LLONG, H5T_IEEE_F64BE, H5T_IEEE_F64LE, H5T_INTEL_B64, - H5T_INTEL_F64, H5T_INTEL_I64, H5T_MIPS_B64, H5T_MIPS_F64, + { H5T_NATIVE_LLONG, H5T_INTEL_B64, + H5T_INTEL_I64, H5T_MIPS_B64, H5T_MIPS_I64, H5T_NATIVE_B64, H5T_NATIVE_INT64, H5T_NATIVE_INT_FAST64, H5T_NATIVE_INT_LEAST64, H5T_STD_B64BE, H5T_STD_B64LE, H5T_STD_I64BE, H5T_STD_I64LE, H5T_UNIX_D64BE, H5T_UNIX_D64LE, H5T_ALPHA_B64, - H5T_ALPHA_F64, H5T_ALPHA_I64 }, + H5T_ALPHA_I64 }, /* GDL_ULONG */ { H5T_NATIVE_ULONG, H5T_ALPHA_U32, H5T_INTEL_U32, H5T_MIPS_U32, @@ -181,9 +181,9 @@ namespace lib { /* GDL_LONG */ { /// H5T_NATIVE_HBOOL, /// ^--- disabled as it matches against 'H5T_STD_U8LE' (GDL_BYTE) - H5T_NATIVE_LONG, H5T_ALPHA_B32, H5T_ALPHA_F32, - H5T_ALPHA_I32, H5T_IEEE_F32BE, H5T_IEEE_F32LE, H5T_INTEL_B32, - H5T_INTEL_F32, H5T_INTEL_I32, H5T_MIPS_B32, H5T_MIPS_F32, + H5T_NATIVE_LONG, H5T_ALPHA_B32, + H5T_ALPHA_I32, H5T_INTEL_B32, + H5T_INTEL_I32, H5T_MIPS_B32, H5T_MIPS_I32, H5T_NATIVE_B32, H5T_NATIVE_INT32, H5T_NATIVE_INT_FAST32, H5T_NATIVE_INT_LEAST32, H5T_STD_B32BE, H5T_STD_B32LE, H5T_STD_I32BE, H5T_STD_I32LE, H5T_UNIX_D32BE, H5T_UNIX_D32LE }, @@ -863,7 +863,8 @@ namespace lib { } else if (ourType == GDL_STRING) { - if (debug) printf("fixed-length string dataset\n"); + bool isVarLenStr = H5Tis_variable_str(elem_dtype) > 0; + if (debug) printf(isVarLenStr ? "variable-length string dataset\n" : "fixed-length string dataset\n"); // string length (terminator included) SizeT str_len = H5Tget_size(elem_dtype); @@ -872,6 +873,17 @@ namespace lib { SizeT num_elems=1; for(int i=0; iGetParGlobal(nParam-1); e->GetParGlobal(nParam-1)=temp; } + else if (att_type == NC_STRING) { +// Read NC_STRING + char **stemp = new char*[length]; +// Read the (possible array of) strings - nc_get_att_string allocates memory for the strings themselves + status = nc_get_att_string(cdfid, varid, attname.c_str(), stemp); + ncdf_handle_error(e, status, "NCDF_ATTGET"); + delete e->GetParGlobal(nParam-1); +// This will only copy first string of array to the returned variable... + e->GetParGlobal(nParam-1)=new DStringGDL(stemp[0]); +// Free the string data from the stemp (having copied the string to the output) + status = nc_free_string(length, stemp); + ncdf_handle_error(e, status, "NCDF_ATTGET"); + } else { dimension dim(length); BaseGDL* temp; switch (att_type) { + case NC_INT64 : + { + long long int *i64p = new long long int[length]; + status=nc_get_att_longlong(cdfid, varid, attname.c_str(), i64p); + ncdf_att_handle_error(e, status, "NCDF_ATTGET", i64p); + temp = length == 1 ? new DLong64GDL(BaseGDL::NOZERO) : new DLong64GDL(dim, BaseGDL::NOZERO); + memcpy(&(*static_cast(temp))[0], &(*i64p), length * sizeof(long long int)); + delete[] i64p; + break; + } + case NC_UINT64 : + { + unsigned long long int *ui64p = new unsigned long long int[length]; + status=nc_get_att_ulonglong(cdfid, varid, attname.c_str(), ui64p); + ncdf_att_handle_error(e, status, "NCDF_ATTGET", ui64p); + temp = length == 1 ? new DULong64GDL(BaseGDL::NOZERO) : new DULong64GDL(dim, BaseGDL::NOZERO); + memcpy(&(*static_cast(temp))[0], &(*ui64p), length * sizeof(unsigned long long int)); + delete[] ui64p; + break; + } case NC_INT : { int *ip = new int[length]; @@ -235,6 +268,16 @@ namespace lib { delete[] ip; break; } + case NC_UINT : + { + unsigned int *uip = new unsigned int[length]; + status=nc_get_att_uint(cdfid, varid, attname.c_str(), uip); + ncdf_att_handle_error(e, status, "NCDF_ATTGET", uip); + temp = length == 1 ? new DULongGDL(BaseGDL::NOZERO) : new DULongGDL(dim, BaseGDL::NOZERO); + memcpy(&(*static_cast(temp))[0], &(*uip), length * sizeof(unsigned int)); + delete[] uip; + break; + } case NC_SHORT : { short *sp = new short[length]; @@ -245,6 +288,16 @@ namespace lib { delete[] sp; break; } + case NC_USHORT : + { + unsigned short *usp = new unsigned short[length]; + status = nc_get_att_ushort(cdfid, varid, attname.c_str(), usp); + ncdf_att_handle_error(e, status, "NCDF_ATTGET", usp); + temp = length == 1 ? new DIntGDL(BaseGDL::NOZERO) : new DIntGDL(dim, BaseGDL::NOZERO); + memcpy(&(*static_cast(temp))[0], &(*usp), length * sizeof(DUInt)); + delete[] usp; + break; + } case NC_FLOAT : { float *fp = new float[length]; @@ -266,9 +319,15 @@ namespace lib { break; } case NC_BYTE : + case NC_UBYTE : { unsigned char *bp = new unsigned char[length]; status = nc_get_att_uchar(cdfid, varid, attname.c_str(), bp); + if (status == NC_ERANGE) { +// just warn on NC_ERANGE (probably intended to if idl char used to represent signed byte) + //Warning("Warning in NCDF_ATTGET: NC_ERANGE while reading BYTE"); + status=NC_NOERR; + } ncdf_att_handle_error(e, status, "NCDF_ATTGET", bp); temp = length == 1 ? new DByteGDL(BaseGDL::NOZERO) : new DByteGDL(dim, BaseGDL::NOZERO); memcpy(&(*static_cast(temp))[0], &(*bp), length * sizeof(DByte)); @@ -330,10 +389,16 @@ namespace lib { if (val->Type() == GDL_BYTE) xtype=NC_BYTE; if (val->Type() == GDL_STRING) xtype=NC_CHAR; if (val->Type() == GDL_INT) xtype=NC_SHORT; + if (val->Type() == GDL_UINT) xtype=NC_USHORT; if (val->Type() == GDL_LONG) xtype=NC_INT; + if (val->Type() == GDL_ULONG) xtype=NC_UINT; + if (val->Type() == GDL_LONG64) xtype=NC_INT64; + if (val->Type() == GDL_ULONG64) xtype=NC_UINT64; if (val->Type() == GDL_FLOAT) xtype=NC_FLOAT; if (val->Type() == GDL_DOUBLE) xtype=NC_DOUBLE; // SA: TODO: GDL_UINT, GDL_ULONG, GDL_COMPLEX, GDL_PTR... + // + // "BYTE","CHAR","DOUBLE","FLOAT","LONG","SHORT","UBYTE","ULONG","USHORT","INT64","UINT64","STRING" if(e->KeywordSet(2)) //GDL_BYTE xtype=NC_BYTE; @@ -347,6 +412,18 @@ namespace lib { xtype=NC_INT; else if(e->KeywordSet(7)) //SHORT xtype=NC_SHORT; + else if(e->KeywordSet(8)) + xtype=NC_UBYTE; + else if(e->KeywordSet(9)) + xtype=NC_UINT; + else if(e->KeywordSet(10)) + xtype=NC_USHORT; + else if(e->KeywordSet(11)) + xtype=NC_INT64; + else if(e->KeywordSet(12)) + xtype=NC_UINT64; + else if(e->KeywordSet(13)) + xtype=NC_STRING; // LENGTH keyword support DLong length; @@ -367,6 +444,11 @@ namespace lib { attname.c_str(),xtype, (size_t)length, (const unsigned char *)&(*bvar)[0]); + if (status == NC_ERANGE) { +// just warn on NC_ERANGE (probably intended to if idl char used to represent signed byte) + //Warning("Warning in NCDF_ATTPUT: NC_ERANGE while writing BYTE"); + status=NC_NOERR; + } } else if(val->Type() == GDL_STRING) { @@ -375,12 +457,17 @@ namespace lib { length = cvar.length(); e->AssureLongScalarKWIfPresent(1, length); - if (length > cvar.length()) e->Throw("LENGTH keyword value (" + i2s(length) + - ") exceedes the data length (" + i2s(cvar.length()) + ")"); - if (length < cvar.length()) cvar.resize(length); - - status=nc_put_att_text(cdfid,varid, attname.c_str(), + if ( xtype == NC_STRING ) { + char *tmp = (char *)cvar.c_str(); + status=nc_put_att_string(cdfid,varid, attname.c_str(), + 1, (const char**)&tmp); + } else { + if (length > cvar.length()) e->Throw("LENGTH keyword value (" + i2s(length) + + ") exceedes the data length (" + i2s(cvar.length()) + ")"); + if (length < cvar.length()) cvar.resize(length); + status=nc_put_att_text(cdfid,varid, attname.c_str(), cvar.length(), (char *)cvar.c_str()); + } } else if(val->Type() == GDL_INT) { @@ -389,6 +476,13 @@ namespace lib { attname.c_str(), xtype, (size_t)length, &(*ivar)[0]); } + else if(val->Type() == GDL_UINT) + { + DUIntGDL * uivar=static_cast(val); + status=nc_put_att_ushort(cdfid,varid, + attname.c_str(), xtype, + (size_t)length, &(*uivar)[0]); + } else if(val->Type() == GDL_LONG) { DLongGDL * lvar=static_cast(val); @@ -396,6 +490,27 @@ namespace lib { attname.c_str(),xtype, (size_t)length, &(*lvar)[0]); } + else if(val->Type() == GDL_ULONG) + { + DULongGDL * ulvar=static_cast(val); + status=nc_put_att_uint(cdfid,varid, + attname.c_str(),xtype, + (size_t)length, &(*ulvar)[0]); + } + else if(val->Type() == GDL_LONG64) + { + DLong64GDL * l64var=static_cast(val); + status=nc_put_att_longlong(cdfid,varid, + attname.c_str(),xtype, + (size_t)length, &(*l64var)[0]); + } + else if(val->Type() == GDL_ULONG64) + { + DULong64GDL * ul64var=static_cast(val); + status=nc_put_att_ulonglong(cdfid,varid, + attname.c_str(),xtype, + (size_t)length, &(*ul64var)[0]); + } else if(val->Type() == GDL_FLOAT) { DFloatGDL * fvar=static_cast(val); @@ -410,7 +525,6 @@ namespace lib { attname.c_str(),xtype, (size_t)length, &(*dvar)[0]); } - ncdf_handle_error(e, status,"NCDF_ATTPUT"); return; diff --git a/src/ncdf_cl.cpp b/src/ncdf_cl.cpp index b8e3c5cef..9caded18b 100644 --- a/src/ncdf_cl.cpp +++ b/src/ncdf_cl.cpp @@ -45,11 +45,15 @@ namespace lib { switch (vartype) { case NC_BYTE: return DStringGDL("BYTE");//8 bit + case NC_UBYTE: return DStringGDL("UBYTE");//8 bit case NC_CHAR: return DStringGDL("CHAR");//8 bit as string case NC_SHORT: return DStringGDL("INT");//16 bit + case NC_USHORT: return DStringGDL("UINT");//16 bit case NC_INT: return DStringGDL("LONG");//32 bit + case NC_UINT: return DStringGDL("ULONG");//32 bit case NC_FLOAT: return DStringGDL("FLOAT");//32 bit float case NC_DOUBLE:return DStringGDL("DOUBLE");//64 bit double + case NC_STRING:return DStringGDL("STRING");//String } return DStringGDL("UNKNOWN"); } diff --git a/src/ncdf_var_cl.cpp b/src/ncdf_var_cl.cpp index cf79e1be7..893eaff9d 100644 --- a/src/ncdf_var_cl.cpp +++ b/src/ncdf_var_cl.cpp @@ -301,6 +301,17 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = new DLongGDL((ivar)); } + else if(var_type == NC_UINT) + { + + unsigned int uivar; + status=nc_get_var1_uint(cdfid, + varid, + index, + &uivar); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = new DULongGDL((uivar)); + } else if(var_type == NC_SHORT) { @@ -312,6 +323,17 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = new DIntGDL((svar)); } + else if(var_type == NC_USHORT) + { + + unsigned short usvar; + status=nc_get_var1_ushort(cdfid, + varid, + index, + &usvar); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = new DUIntGDL((usvar)); + } else if(var_type == NC_CHAR) { char cvar; @@ -319,7 +341,7 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = new DByteGDL(cvar); } - else if(var_type == NC_BYTE){ + else if(var_type == NC_BYTE || var_type == NC_UBYTE ){ unsigned char bvar; status=nc_get_var1_uchar(cdfid, varid, @@ -436,6 +458,14 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2)=temp; } + else if (var_type == NC_USHORT) + { + DUIntGDL * temp=new DUIntGDL(dim,BaseGDL::NOZERO); + status=nc_get_var_ushort(cdfid, varid,&(*temp)[0]); + ncdf_var_handle_error(e, status, "NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2)=temp; + } else if (var_type == NC_INT) { DLongGDL* temp=new DLongGDL(dim,BaseGDL::NOZERO); @@ -444,15 +474,26 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2)=temp; } - else if (var_type == NC_BYTE) + else if (var_type == NC_UINT) + { + DULongGDL* temp=new DULongGDL(dim,BaseGDL::NOZERO); + status=nc_get_var_uint(cdfid, varid,&(*temp)[0]); + ncdf_var_handle_error(e, status, "NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2)=temp; + } + else if (var_type == NC_BYTE || var_type == NC_UBYTE ) { DByteGDL* temp=new DByteGDL(dim,BaseGDL::NOZERO); status=nc_get_var_uchar(cdfid, varid, &(*temp)[0]); if (status != NC_ERANGE) { ncdf_var_handle_error(e,status,"NCDF_VARGET (ici)", temp); } else { - Warning("Warning in NCDF_VARGET: NC_ERANGE during BYTE reading"); - ncdf_var_handle_error(e,status,"NCDF_VARGET (ici)", temp); +// exceeding range of signed byte is allowed in IDL +// because IDL has no signed byte datatype - need to read it as unsigned byte and +// then code needs to handle the sign bit (convert to int then values>127 need 256 subtracting) + //Warning("Warning in NCDF_VARGET: NC_ERANGE while reading BYTE"); + status=NC_NOERR; } GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2)=temp; @@ -563,6 +604,14 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = temp; } + else if (var_type == NC_USHORT) + { + DUIntGDL *temp = new DUIntGDL(dim, BaseGDL::NOZERO); + status = nc_get_vara_ushort(cdfid, varid, off, cou, &(*temp)[0]); + ncdf_var_handle_error(e, status, "NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = temp; + } else if(var_type == NC_INT) { DLongGDL *temp = new DLongGDL(dim,BaseGDL::NOZERO); @@ -571,7 +620,15 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = temp; } - else if (var_type == NC_BYTE) + else if(var_type == NC_UINT) + { + DULongGDL *temp = new DULongGDL(dim,BaseGDL::NOZERO); + status = nc_get_vara_uint(cdfid, varid, off, cou, &(*temp)[0]); + ncdf_var_handle_error(e, status, "NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = temp; + } + else if (var_type == NC_BYTE || var_type == NC_UBYTE ) { DByteGDL *temp=new DByteGDL(dim,BaseGDL::NOZERO); status = nc_get_vara_uchar(cdfid, varid, off, cou, &(*temp)[0]); @@ -666,6 +723,14 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = temp; } + else if (var_type == NC_USHORT) + { + DUIntGDL *temp = new DUIntGDL(dim, BaseGDL::NOZERO); + status = nc_get_vars_ushort(cdfid, varid, off, cou, stri, &(*temp)[0]); + ncdf_var_handle_error(e,status,"NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = temp; + } else if(var_type == NC_INT) { DLongGDL *temp = new DLongGDL(dim, BaseGDL::NOZERO); @@ -674,7 +739,15 @@ else if(var_type == NC_LONG) GDLDelete(e->GetParGlobal(2)); e->GetParGlobal(2) = temp; } - else if(var_type == NC_BYTE) + else if(var_type == NC_UINT) + { + DULongGDL *temp = new DULongGDL(dim, BaseGDL::NOZERO); + status = nc_get_vars_uint(cdfid, varid, off,cou, stri, &(*temp)[0]); + ncdf_var_handle_error(e, status, "NCDF_VARGET", temp); + GDLDelete(e->GetParGlobal(2)); + e->GetParGlobal(2) = temp; + } + else if(var_type == NC_BYTE || var_type == NC_UBYTE ) { DByteGDL *temp=new DByteGDL(dim, BaseGDL::NOZERO); status = nc_get_vars_uchar(cdfid, varid, off, cou, stri, &(*temp)[0]); @@ -734,6 +807,7 @@ else if(var_type == NC_LONG) //dims is not set, scalar } +//"BYTE","CHAR","DOUBLE","FLOAT","LONG","SHORT","UBYTE","ULONG","USHORT","INT64","UINT64","STRING","GZIP","SHUFFLE" if(e->KeywordSet(0))//GDL_BYTE type=NC_BYTE; else if(e->KeywordSet(1))//CHAR @@ -744,6 +818,18 @@ else if(var_type == NC_LONG) type=NC_INT; else if(e->KeywordSet(5))//SHORT type=NC_SHORT; + else if(e->KeywordSet(6)) + type=NC_UBYTE; + else if(e->KeywordSet(7)) + type=NC_UINT; + else if(e->KeywordSet(8)) + type=NC_USHORT; + else if(e->KeywordSet(9)) + type=NC_INT64; + else if(e->KeywordSet(10)) + type=NC_UINT64; + else if(e->KeywordSet(11)) + type=NC_STRING; else type=NC_FLOAT; @@ -759,6 +845,20 @@ else if(var_type == NC_LONG) "Unable to define variable, name in use by another variable ("+var_name+")"); else ncdf_handle_error(e,status,"NCDF_VARDEF"); + +// Compression settings (controlled by GZIP keyword) + int deflate_level = 0; + int shuffle = NC_NOSHUFFLE; + if(e->KeywordSet(12)) { + DLong deflLong = 0; + e->AssureLongScalarKW(12, deflLong); + deflate_level = (int) deflLong; + if(e->KeywordSet(13)) shuffle=NC_SHUFFLE; +//printf("Compression on level=%d; shuffle=%d\n",deflate_level,shuffle); + status = nc_def_var_deflate(cdfid, var_id, shuffle, 1, deflate_level); + ncdf_handle_error(e,status,"NCDF_VARDEF"); + } + return new DIntGDL(var_id); @@ -1031,9 +1131,23 @@ else if(var_type == NC_LONG) break; } case GDL_STRING : - status = nc_put_vars_text(cdfid, varid, offset, count, stride, - (*static_cast(v))[0].c_str()); + { + const char *tmp = (*static_cast(v))[0].c_str(); + switch (var_type) + { + case NC_CHAR : // char array + status = nc_put_vars_text(cdfid, varid, offset, count, stride, tmp); + break; + case NC_STRING : // string +//printf("NC_STRING %d %d\n",var_ndims,count[0]); + if ( var_ndims > 1 || count[0] != 1) e->Throw("GDL string array cannot be written - only single value allowed by code!"); + status = nc_put_vars_string(cdfid, varid, offset, count, stride, &tmp); + break; + default: + e->Throw("GDL string can only be written to NC_CHAR or NC_STRING!"); + } break; + } // reporting illegal types (could be done before...) case GDL_STRUCT : e->Throw("Struct expression not allowed in this context: " diff --git a/src/pro/map/map_proj_init.pro b/src/pro/map/map_proj_init.pro index 5636fa07e..9af79f976 100644 --- a/src/pro/map/map_proj_init.pro +++ b/src/pro/map/map_proj_init.pro @@ -577,7 +577,8 @@ endif else begin ; not interrupted ;alpha is the rotation angle for ;oblique cylindric. alpha=0 for ;transverse projections - val=80 +; RS fix for standard map projection which should not be cut at 80deg lat + if index eq 27 then val=90 else val=80 ; remove poles map_clip_set, map=myMap, clip_plane=[myMap.pole[4:6],sin(deg2rad*val)] map_clip_set, map=myMap, clip_plane=[-1*myMap.pole[4:6],sin(deg2rad*val)]