[ncl-talk] Visualization of SAF NWC products
Rick Brownrigg
brownrig at ucar.edu
Mon Feb 22 10:04:07 MST 2016
Hi,
There are a couple of issues:
i) I think perhaps you have rows/columns transposed in calculating your x,y
cells. I I replace the two do-loops with the following, I get coordinates
that span the range reported in the file's metadata:
do i=0,numCols-1
x(:,i) = projectedMinX + i*cellSize
end do
do i=0,numRows-1
y(i,:) = projectedMinY + i*cellSize
end do
ii) At the maximum extremes of the cells, proj4 reports back "infinity" --
it perhaps represents a discontinuity in this particular projection? At
any rate, if I replace thos INFs with _FillValue, I then get a map that
looks like your previous map over Greece. Right after the call to
transform_coordinate, try adding:
x = where(x.gt.1000., x at _FillValue, x)
y = where(y.gt.1000., y at _FillValue, y)
Hope that helps...
Rick
On Mon, Feb 22, 2016 at 4:14 AM, Athanasios Karagiannidis <thankar at live.com>
wrote:
> Dear Rick
>
> I encountered a new problem when trying to adapt the code for a bigger
> domain.
>
> To begin with, when using the code for the Greek area image, everything
> works just fine.
>
> I then tried to use the code for the visualization of a bigger domain.
> The code could did not compute correctly the maximum latitude and longitude
> and a chain of errors followed. Please see below:
>
> [safnwc at 192 example_scripts]$ ncl geos2LL_europe.ncl
> Copyright (C) 1995-2014 - All Rights Reserved
> University Corporation for Atmospheric Research
> NCAR Command Language Version 6.2.1
> The use of this software is governed by a License Agreement.
> See http://www.ncl.ucar.edu/ for more details.
> (0) Span in X direction (meters): 3.59748e+06
> (0) Span in Y direction (meters): 2.33731e+06
> (0) calculated maxX: 2.77237e+06 vs. reported in file: 2.77237e+06
> (0) calculated maxY: 5.05868e+06 vs. reported in file: 5.05868e+06
> (0) longitude range: -42.10171820490142, inf
> (0) latitude range: 25.96460103290488, inf
> (0) gsn_csm_map_ce: Warning: you set mpMaxLonF to a value > 180, but
> (0) didn't set mpCenterLonF. Setting mpCenterLonF to inf
> warning:tofloat: there are 1 double larger than FLT_MAX, which has been
> flagged missing.
> warning:tofloat: there are 1 double larger than FLT_MAX, which has been
> flagged missing.
> fatal:Subscript out of range, error in subscript #0
> fatal:An error occurred reading wharray
> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 5843 in
> file /home/safnwc/ncl-6.2.1/lib/ncarg/nclscripts/csm/gsn_csm.ncl
>
> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 6639 in
> file /home/safnwc/ncl-6.2.1/lib/ncarg/nclscripts/csm/gsn_csm.ncl
>
> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 61 in
> file geos2LL_europe.ncl
>
> fatal:Variable (map) is undefined
> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 70 in
> file geos2LL_europe.ncl
>
> fatal:Variable (map) is undefined
> fatal:["Execute.c":8578]:Execute: Error occurred at or near line 72 in
> file geos2LL_europe.nc
>
>
> I send you the hdf file of the new domain, and my adapted code. I changed
> the definitions of x, y and z matrices to (/ Lines Cols/) from (/Cols
> Lines/) in your original script, and also the associated do loop limits
> just below them. Of course, i also used the European domain parameters.
>
> I fail to see the reason for that failure, since the code for the smaller
> domain work perfectly.
>
> Kind regards
>
> Sakis
>
>
> Athanasios F. Karagiannidis
> Physicist, PhD Meteorology and Climatology
> E-mail: thankar at live.com
>
>
>
> ------------------------------
> Date: Wed, 17 Feb 2016 07:19:14 -0700
> Subject: Re: [ncl-talk] Visualization of SAF NWC products
> From: brownrig at ucar.edu
> To: thankar at live.com
> CC: ncl-talk at ucar.edu
>
>
> 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/20160222/0c80ccef/attachment.html
More information about the ncl-talk
mailing list