[ncl-talk] Visualization of SAF NWC products
Rick Brownrigg
brownrig at ucar.edu
Wed Feb 17 07:19:14 MST 2016
Hi Sakis,
I have determined it is indeed a proj4 version issue, but fortunately with
an easy to deploy fix. In stepping through the code, I discovered that the
ellipse *must* be specified for lonlat projections, instead of defaulting
to spherical as in previous versions of proj4 (although I can find nothing
in the docs or ChangeLog that spells this out).
To fix in your script, replace the 2nd line that reads:
dstProj = "+proj=lonlat"
with
dstProj = "+proj=lonlat +ellps=sphere"
That should do it...
Rick
On Tue, Feb 16, 2016 at 10:46 PM, Rick Brownrigg <brownrig at ucar.edu> wrote:
> Hi,
>
> I have been able to verify on my end that this problem is version specific
> to proj4. The script works fine on the slightly older 4.8.0 versions, but
> appears in the latest 4.9.1 proj4 (Sakis, can you verify this on your side
> -- what version gets reported if you simply type "proj" on the command
> line?). This seems like a bug in the proj4 library, but I can find no
> mention of it in their issues database. I will try to gather more
> information and file a bug report with the proj4 folks.
>
> In the meantime, Sakis, are you able to down-version your proj4 library to
> something like 4.8.0?
>
> Rick
>
> On Mon, Feb 15, 2016 at 2:46 AM, Athanasios Karagiannidis <
> thankar at live.com> wrote:
>
>> Dear Rick,
>>
>> I understand that you did not manage to identify the problem.
>> If that is the problem , are there any suggestions regarding the
>> visualization of SAF NWC (basically meteosat) products?
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> Date: Mon, 8 Feb 2016 20:38:47 -0700
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> From: brownrig at ucar.edu
>> To: thankar at live.com
>>
>>
>> Hi Sakis,
>>
>> Offhand I don't know what is going on here. It almost appears that the
>> proverbial "rug has been yanked out from under us", with respect to the
>> projection software -- I've used that package for over a dozen years, and
>> "+proj=latlon" is/was by all accounts a legitimate projection spec. I get
>> the same result as you on my system, which has the latest proj4 release --
>> I want to try and track down a system with an older release in the morning
>> and see if I get the same results.
>>
>> Rick
>>
>> On Mon, Feb 8, 2016 at 9:10 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Dear Rick
>>
>> I'm coming back with a problem regarding an earlier problem i had some
>> time ago.
>>
>> I was trying to visualize Meteosat projection data (for example the
>> CTTH_HEIGHT parameter in the h5 file) and you provided the script i attach
>> which is supposed to draw the png i also attach.
>>
>> However when i try to execute, i get the following message:
>>
>> (0) Span in X direction (meters): 2.16929e+06
>> (0) Span in Y direction (meters): 1.53321e+06
>> (0) calculated maxX: 3.03041e+06 vs. reported in file: 3.03041e+06
>> (0) calculated maxY: 4.4706e+06 vs. reported in file: 4.4706e+06
>> FATAL ERROR: Can not make a projection out of <+proj=lonlat>
>> fatal:TransformCoordinate: Couldn't transform coordinates.
>> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 38 in
>> file geos2LL.ncl
>>
>> (0) longitude range: 861115.75, 3030407.25
>> (0) latitude range: 2937395, 4470601
>> (0) gsn_csm_map_ce: Warning: you set mpMaxLonF to a value > 180, but
>> (0) didn't set mpCenterLonF. Setting mpCenterLonF to
>> 1945761.5
>> warning:MapTransInitialize: mpMinLonF out of range: resetting
>> warning:MapTransInitialize: mpMaxLonF out of range: resetting
>> warning:MapSetTrans: map limits invalid - using maximal area
>> warning:mpMapLimitMode is not a valid resource in map at this time
>> (0) get_lat_values: Warning: Your latitude values do not fall between
>> -90 and 90 inclusive.
>> (0) You will not get 'nice' latitude labels.
>> warning:Attempt to reference attribute (minor) which is undefined
>>
>> I think the problem begins at *dstProj = "+proj=lonlat"*.
>> As a result nothing is plotted.
>> Could you help me?
>>
>> I attach a data file in case you need it.
>>
>> Thank you in advance
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> Date: Fri, 5 Jun 2015 15:16:34 -0600
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> From: brownrig at ucar.edu
>> To: thankar at live.com
>> CC: haley at ucar.edu
>>
>> Hi Sakis,
>>
>> The file does indeed seem to have all the requisite projection
>> information in it to calculate latlon coordinates, using, our interface to
>> the proj4 library. I've attached a script that does just that. I've also
>> attached the the plot it generates. I'm a bit troubled by the shape of the
>> lat/lon grid, and by the fact that the spans in X/Y in meters do not align
>> with what gets reported in Google Earth. However, the native projection is
>> "geostationary view", and its possible that it imparts rather severe
>> distortion of distance at the edges of the projection, so maybe the values
>> are reasonable.
>>
>> Do you have a feel for what the shape/extent of your data should look
>> like?
>>
>> Let me know what you think...
>> Rick
>>
>>
>>
>> On Fri, Jun 5, 2015 at 12:53 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Dear Mary
>>
>> Thank you for your help!
>>
>>
>> The speed of the code was indeed an issue.
>>
>> Regarding the lat / lon issue:
>> 1. Yes, they are not equally spaced.
>> 2. Yes, the whole information about the way the rasterized data
>> should be assigned to lat/lon is (to my knowledge) included in the h5 files
>> attributes.
>> 3. I have a way to bypass the problem by loading 2 separate files
>> that provide lat/lon info. However I would like to be independent of them,
>> for a series of reasons.
>> I’m looking forward to Rick’s response.
>>
>> Kind regards
>>
>> Sakis
>>
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> Date: Thu, 4 Jun 2015 22:59:55 -0600
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> From: haley at ucar.edu
>> To: thankar at live.com
>> CC: brownrig at ucar.edu
>>
>>
>> Dear Sakis,
>>
>> You've been so patient, thanks! As a start, I recommend that you use
>> raster mode, as it is much faster:
>>
>> res at cnFillMode = "RasterFill"
>>
>> The code ran in about 1 second, compared to minutes before.
>>
>> The issue here is that you don't have lat/lon arrays on the file, so you
>> have to either construct them yourself, as you are trying to do with the
>> fspan function, or else if you know the exact map projection, then you can
>> just draw the contours in the correct map project, with res at tfDoNDCOverlay
>> = True. This is basically telling NCL that you know what the projection
>> is, and to not require lat/lon coordinates.
>>
>> I think the use of "fspan" to generate the lat/lon values may be too
>> simplistic, as this is assuming your lat/lon values are equally spaced, and
>> not rotated. I have a feeling that neither of these conditions are true
>> for your data. But, I could be wrong! I just know that the plot doesn't
>> look quite right if you run the script as it.
>>
>> I can't tell from your file what the correct map projection is, as I'm
>> not familiar with the terminology it's using. It looks like the information
>> is all there in the file attribute information:
>> GEOTRANSFORM_GDAL_TABLE : -5570248.832537, 3000.403357, 0.000000,
>> 5570248.832537,0.000000, -3000.403357
>> PROJECTION : +proj=geos +a=6378169.0 +b=6356583.8 +lon_0=0.0
>> +h=35785831.000, -3000.403357
>>
>> PROJECTION_NAME : GEOS<+000.0>a=6378169.0 +b=6356583.8 +lon_0=0.0
>> +h=35785831.000, -3000.403357
>> REGION_NAME : Greece000.0>a=6378169.0 +b=6356583.8 +lon_0=0.0
>> +h=35785831.000, -3000.403357
>> XGEO_LOW_RIGHT : 3030407.390714755
>>
>> XGEO_UP_LEFT : 861115.7635001335
>> YGEO_LOW_RIGHT : 2937394.886643311
>> YGEO_UP_LEFT : 4470601.002143549
>>
>> My guess is that you might be able to use one of our unadvertised,
>> internal GDAL transformation routines to convert this information into
>> lat/lon coordinates. I've CC'ed Rick Brownrigg who may know something
>> about this.
>>
>>
>> --Mary
>>
>>
>>
>> On Sun, May 31, 2015 at 1:48 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Dear Mary
>>
>> I attached the latets script that was prepared for me along with a sample
>> dataset (h5 file).
>>
>> Thank you for your interest and help.
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> Date: Sat, 30 May 2015 15:34:05 -0600
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> From: haley at ucar.edu
>> To: thankar at live.com
>> CC: ncl-talk at ucar.edu
>>
>>
>> Dear Sakis,
>>
>> Can you provide the data and the latest script that you are using, so we
>> can take a look? You can use our ftp if you don't have a place to put it:
>>
>> http://www.ncl.ucar.edu/report_bug.shtml#HowToFTP
>>
>> --Mary
>>
>>
>> On Thu, May 28, 2015 at 2:05 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Good morning to all!
>>
>> Following the earlier conversation (to which I respond) I think I have
>> the obligation to clarify some issues, since these conversations are
>> available to all users, and they are used to solve various kinds of
>> problems.
>> After doing some more tests, I found out that the otherwise elegant and
>> technically flawless script that was kindly prepared for me, does not
>> illustrates correctly the CTTH_HEIGHT field. The problem is related to the
>> fact that the line/column space of the raster data does not refer to a
>> regular lat/lon space. In simple words the pixels of the raster data are
>> not equal in lat/lon dimensions.
>>
>> Regarding more technical issues:
>> XGEO_LOW_RIGHT, XGEO_UP_LEFT, YGEO_LOW_RIGHT, YGEO_UP_LEFT are not given
>> in degrees but in meters.
>>
>> dims = dimsizes(CT) , lat = fspan(minLat, maxLat, dims(0)) and lon =
>> fspan(minLon, maxLon, dims(1)) are producing a regular grid (which is
>> not correct)
>>
>> It would be really helpful if a script could project correctly Meteosat
>> fields (I think safnwc is using the Meteosat grid) but I’m not sure that
>> such a thing exists.
>>
>> Kind regards
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> From: thankar at live.com
>> To: haley at ucar.edu
>> Date: Sat, 21 Feb 2015 08:29:02 +0000
>> CC: ncl-talk at ucar.edu
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>>
>> Mary, thank you for your input!
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> Date: Fri, 20 Feb 2015 09:46:02 -0700
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> From: haley at ucar.edu
>> To: thankar at live.com
>> CC: huangwei at ucar.edu; ncl-talk at ucar.edu
>>
>> The issue with "namedcolor2rgb" is not really an issue. The code should
>> still work. The error message is just misleading, but can be ignored.
>>
>> We have fixed the incorrect error message for the V6.3.0 release coming
>> next month.
>>
>> --Mary
>>
>>
>> On Fri, Feb 20, 2015 at 1:46 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Wei,
>>
>> Your help is more than appreciated.
>>
>> I will update my NCL version.
>>
>> Kind regards
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> From: huangwei at ucar.edu
>> Date: Thu, 19 Feb 2015 09:23:14 -0700
>>
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> To: thankar at live.com
>> CC: ncl-talk at ucar.edu
>>
>> Sakis,
>>
>> Attached is script which looks more like what you see with HDFview,
>> but used a whole user-defined color table.
>>
>> "gsdtol" is a standard color-table name, in NCL 6.2.1. Your NCL of 6.1.2
>> is a little bit old, and I strongly suggest you update your NCL to the
>> newest version.
>>
>> Compare the newly attached script and the script before,
>> you'll be able to figure out that the difference is the projection
>> applied.
>> So you can turn on/off those options based on your requirement.
>>
>> A new png file is attached as well.
>>
>> Regards,
>>
>> Wei
>>
>>
>>
>> ================================================
>> 1850 Table Mesa Dr.
>> Boulder, CO 80307
>> Phone: 303-497-8924
>>
>> On Thu, Feb 19, 2015 at 4:07 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Dear Wei
>>
>> I'm more than thankful for your script! It is excellent in correctly
>> reading my data and more or less what i wanted regarding the visualization.
>> I have just 3 comments:
>> 1. I would really like to see the drawn field (in example the CT) even
>> over land. How can i do this?
>> 2. When i run the script in my machine i get the following warnings:
>>
>>
>>
>> *(0) Warning: namedcolor2rgb: 'gsdtol' is not a valid named
>> color.(0) Will return missing values for this color.*
>> Is 'gsdtol' a custom made color? In this case maybe i should replace it.
>>
>> 3. In the example HDFView visualization i provided the grid was not
>> "tilted" to the left (a projection issue i guess!). How should i make your
>> visualization, so it looks like the HDFView one?
>>
>> Thank you in advance!
>>
>> Sakis
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>>
>> ------------------------------
>> From: huangwei at ucar.edu
>> Date: Wed, 18 Feb 2015 13:26:21 -0700
>> Subject: Re: [ncl-talk] Visualization of SAF NWC products
>> To: thankar at live.com
>> CC: ncl-talk at ucar.edu
>>
>>
>> Sakis,
>>
>> You may take a look of the attached script and pdf file to see if these
>> are what you wanted.
>>
>> Thanks,
>>
>> Wei
>>
>>
>>
>> ================================================
>> 1850 Table Mesa Dr.
>> Boulder, CO 80307
>> Phone: 303-497-8924
>>
>> On Wed, Feb 18, 2015 at 1:38 AM, Athanasios Karagiannidis <
>> thankar at live.com> wrote:
>>
>> Dear sirs/madams
>>
>> I want to visualize SAF NWC output data using ncl. They are included in
>> an h5 file.
>> These data are not in a simple lat-lon grid.
>> After an earlier contact with you i managed to visualize them by creating
>> such a grid, but i feel that the image is not accurate.
>> What i want to help me with, is to tell me weather there is a standard
>> way or an already existing script to visualize correctly that kind of data.
>> To help you I provide the following:
>>
>> 1. A sample datafile.
>> 2. The printVarSummary output
>> 3. A screenshot of the HDFView visualization of the variable where the
>> Variable info are shown (more or less the same info with printVarSummary).
>> 4. The script i used (although, i do not think that you need that!). The
>> lat/lon limits that i have n my script are not acurate. I used them to make
>> the visualization.
>>
>>
>> ---------- printVarSummary------------------------------
>>
>> ncl 2> printVarSummary(CT)
>>
>> Variable: CT
>> Type: ubyte
>> Total Size: 370688 bytes
>> 370688 values
>> Number of Dimensions: 2
>> Dimensions and sizes: [DIM_000 | 512] x [DIM_001 | 724]
>> Coordinates:
>> Number Of Attributes: 11
>> CLASS : IMAGE_INDEXED
>> ID : CTTH_HEIGHTED
>> IMAGE_COLORMODEL : RGBH_HEIGHTED
>> IMAGE_SUBCLASS : IMAGE_INDEXED
>> IMAGE_VERSION : 1.0GE_INDEXED
>> N_COLS : 724
>> N_LINES : 512
>> OFFSET : -2000
>> PALETTE : ��8�3
>> PRODUCT : CTTHE_INDEXED
>> SCALING_FACTOR : 200
>>
>>
>> ---------------------------------------------------------------------------------
>>
>> ---------------- My Script
>> ---------------------------------------------------
>>
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>
>> begin
>> cdf_file = addfile("SAFNWC_MSG2_CTTH_200812121100_Greece______.h5","r")
>> CT = cdf_file->CTTH_HEIGHT(:,:)
>>
>> res = True
>>
>> lat = fspan(50,20,512) ; A.K. -- reverse the matrix
>> lon = fspan(10,30,724)
>> CT!0 = "lat" ;-- set name of CT dimension 0
>> CT!1 = "lon" ;-- set name of CT dimension 1
>> CT&lat = lat ;-- assign lat array to CT
>> CT&lon = lon ;-- assign lon array to CT
>> CT&lat at units = "degrees_north" ;-- set the units of the named
>> coordinate lat
>> CT&lon at units = "degrees_east" ;-- set the units of the named
>> coordinate lon
>>
>> res at gsnAddCyclic = False
>> res at tiMainString = "Default Color" ; main title
>> res at cnFillOn = True ; turn on color fill
>> res at mpMaxLatF = 50 ; choose
>> subregion
>> res at mpMinLatF = 20
>> res at mpMaxLonF = 30
>> res at mpMinLonF = 10
>>
>> xwks = gsn_open_wks("pdf","CloudTops")
>> plot = gsn_csm_contour_map(xwks,CT,res)
>> end
>>
>>
>> --------------------------------------------------------------------------------
>>
>> Thank you in advance.
>>
>>
>> Sakis
>>
>>
>>
>>
>> Athanasios F. Karagiannidis
>> Physicist, PhD Meteorology and Climatology
>> E-mail: thankar at live.com
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>>
>> _______________________________________________ ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
>>
>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160217/10d0caa9/attachment.html
More information about the ncl-talk
mailing list