Re: [mapserver-users] GRIB files - gdal_translate

2009-05-21 Thread Roberto Garcia,MSc

Hi Frank, tks for the answer.

CTL is an ASCII file that follows Grib files (at least here) and informs 
the metadata of the Grib file. If anyone consider  helping me  my data 
is available by ftp in ftp1.cptec.inpe.br inside bdados folder, as my 
map file is.


I access them with
http://localhost/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/grib.mapmode=maplayer=grib1
or
http://localhost/cgi-bin/mapserv.exe?map=/ms4w/apps/tutorial/htdocs/grib.mapmode=maplayer=tif1

And the results are almost the same (the diff is the color), a line 
below S.A.


Tks

Frank Warmerdam wrote:

Roberto Garcia,MSc wrote:

Size is *161, 98*
Coordinate System is:
.
Origin = (-150.000,-61.2460002)
Pixel Size = (0.938,-0.001836734693878)

...

Band 1 Block=*161x1* Type=Float64, ColorInterp=Undefined


As u can see the size according to the its header is 161x98, but all 
the bands are related as 161x1 block. Why does it happen? 


Roberto,

This means that the file (and bands) are 161x98, but the bands
are subdivided into blocks which are 161x1 implying that the natural
access chunk is a scanline.  This is normal.

Another thing I noticed is that the dif between upper and lower 
latitude in the corner coordinates is very small (-61.246 to -61.426).


This appears to be screwed up.  Perhaps you could provide the file,
and file a ticket about this?

The issue is that, when mapserver renders this layer it only draws a 
straight line below South America, exactly where the whole image 
should start, as my CTL below shows.


What is CTL?

PS. this would likely be better addressed on gdal-dev.  Any ticket
should be filed in the GDAL Trac.

  http://trac.osgeo.org/gdal/

Best regards,


dset ^T126L28.grb
title pac/sa/atlan sector: troposphere file
undef 1e+20
dtype grib
index ^T126L28.gmp
xdef 161 linear -150.15 0.937500
ydef 98 levels
-61.246 -60.311
-59.376 -58.441 -57.506 -56.571 -55.636 -54.701 -53.766 -52.831 
-51.896 -50.961
-50.026 -49.091 -48.156 -47.221 -46.286 -45.350 -44.415 -43.480 
-42.545 -41.610
-40.675 -39.740 -38.805 -37.870 -36.935 -36.000 -35.065 -34.130 
-33.195 -32.260
-31.325 -30.389 -29.454 -28.519 -27.584 -26.649 -25.714 -24.779 
-23.844 -22.909
-21.974 -21.039 -20.104 -19.169 -18.234 -17.299 -16.364 -15.429 
-14.493 -13.558
-12.623 -11.688 -10.753  -9.818  -8.883  -7.948  -7.013  -6.078  
-5.143  -4.208
 -3.273  -2.338  -1.403  -0.468   0.468   1.403   2.338   3.273   
4.208   5.143
  6.078   7.013   7.948   8.883   9.818  10.753  11.688  12.623  
13.558  14.493
 15.429  16.364  17.299  18.234  19.169  20.104  21.039  21.974  
22.909  23.844

 24.779  25.714  26.649  27.584  28.519  29.454
zdef 7 levels
1000 925 850 700 500 300 200   === set Z 1 ou 2 ou 3 ... ou 7 (GRADS)
tdef 1 linear 6Z16may2009  6hr
vars 15
psnm  02,102,  0,  0 SEA LEVEL PRESSURE [hPa]
tsfc  0  187,  1,  0,  0 SURFACE TEMPERATURE [K]
prec  0   61,  1,  0,  0 TOTAL PRECIPITATION [Kg/m2/day]
cbnv  0   71,  3,  0,  0 CLOUD COVER [0-1]
role  0  114,  8,  0,  0 OUTGOING LONG WAVE AT TOP [W/m2]
mask  7  137,100 MASK [-/+]
uvel  7   33,100 ZONAL WIND (U) [m/s]
vvel  7   34,100 MERIDIONAL WIND (V) [m/s]
omeg  7   39,100 OMEGA [Pa/s]
fcor  7   35,100 STREAM FUNCTION [m2/s]
potv  7   36,100 VELOCITY POTENTIAL [m2/s]
zgeo  77,100 GEOPOTENTIAL HEIGHT [gpm]
temp  7   11,100 ABSOLUTE TEMPERATURE [K]
umrl  7   52,100 RELATIVE HUMIDITY [no Dim]
umes  7   51,100 SPECIFIC HUMIDITY [kg/kg]
endvars

Can someone help me to find out what is happening?

Tks a lot

Roberto Garcia


Frank Warmerdam wrote:

Roberto Garcia,MSc wrote:

Hi, tks for the answers.

According to my gdainfo output file attached, can anyone tell me 
what element could be a CLASSITEM?


Roberto,

In MapServer you always classify rasters based on the item [pixel] but
you do not need to explicitly state a classitem.

This is a relatively simple example of raster classification.

LAYER
  NAME grid1
  TYPE raster
  STATUS default
  DATA data/float.tif
  CLASS
NAME red
EXPRESSION ([pixel]  -3)
COLOR 255 0 0
  END
  CLASS
NAME green
EXPRESSION ([pixel] = -3 and [pixel]  3)
COLOR 0 255 0
  END
  CLASS
NAME blue
EXPRESSION ([pixel] = 3)
COLOR 0 0 255
  END
END

I can see my grib using Grads but what tool do I use to translate 
individual bands into GeoTIFF?


You can use gdal_translate to extract particular bands from a grib file
as a geotiff (or all the bands).  For instance, the following would
extract band 3 (total precipitation):

gdal_translate -b 3 T126L28.grb total_precip.tif

Converting the data on-the-fly to an image is the right way to work? 


Having MapServer access the grib file on the fly rather than depending
on pre-translated extracts will reduce data duplication and may help
avoid confusion about what is up to date.  However, it is possible that
there will be a performance penalty for extracting directly from GRIB.
If so, it may

[mapserver-users] GRIB files - gdal_translate

2009-05-20 Thread Roberto Garcia,MSc

Hi people,

I used gdal_translate in the GRIB file I want to show with mapserver, 
but there still is something wrong. The following is a partial output 
from it:


Driver: GRIB/GRIdded Binary (.grb)
Files: data\raster\T126L28.grb
Size is *161, 98*
Coordinate System is:
.
Origin = (-150.000,-61.2460002)
Pixel Size = (0.938,-0.001836734693878)
Corner Coordinates:
Upper Left  (-150.000, -61.246) (150d 0'0.00W, 61d14'45.60S)
Lower Left  (-150.000, -61.426) (150d 0'0.00W, 61d25'33.60S)
Upper Right (   1.018, -61.246) (  1d 1'4.80E, 61d14'45.60S)
Lower Right (   1.018, -61.426) (  1d 1'4.80E, 61d25'33.60S)
Center  ( -74.491, -61.336) ( 74d29'27.60W, 61d20'9.60S)
Band 1 Block=*161x1* Type=Float64, ColorInterp=Undefined
 Description = 0[-] MSL (Mean sea level)
 Metadata:
   GRIB_UNIT=[Pa]
   GRIB_COMMENT=Pressure reduced to MSL [Pa]
   GRIB_ELEMENT=PRMSL
   GRIB_SHORT_NAME=0-MSL
   GRIB_REF_TIME=  1242453600 sec UTC
   GRIB_VALID_TIME=  1242453600 sec UTC
   GRIB_FORECAST_SECONDS=0 sec


As u can see the size according to the its header is 161x98, but all the 
bands are related as 161x1 block. Why does it happen? Another thing I 
noticed is that the dif between upper and lower latitude in the corner 
coordinates is very small (-61.246 to -61.426).


The issue is that, when mapserver renders this layer it only draws a 
straight line below South America, exactly where the whole image should 
start, as my CTL below shows.


dset ^T126L28.grb
title pac/sa/atlan sector: troposphere file
undef 1e+20
dtype grib
index ^T126L28.gmp
xdef 161 linear -150.15 0.937500
ydef 98 levels
-61.246 -60.311
-59.376 -58.441 -57.506 -56.571 -55.636 -54.701 -53.766 -52.831 -51.896 
-50.961
-50.026 -49.091 -48.156 -47.221 -46.286 -45.350 -44.415 -43.480 -42.545 
-41.610
-40.675 -39.740 -38.805 -37.870 -36.935 -36.000 -35.065 -34.130 -33.195 
-32.260
-31.325 -30.389 -29.454 -28.519 -27.584 -26.649 -25.714 -24.779 -23.844 
-22.909
-21.974 -21.039 -20.104 -19.169 -18.234 -17.299 -16.364 -15.429 -14.493 
-13.558
-12.623 -11.688 -10.753  -9.818  -8.883  -7.948  -7.013  -6.078  -5.143  
-4.208
-3.273  -2.338  -1.403  -0.468   0.468   1.403   2.338   3.273   
4.208   5.143
 6.078   7.013   7.948   8.883   9.818  10.753  11.688  12.623  13.558  
14.493
15.429  16.364  17.299  18.234  19.169  20.104  21.039  21.974  22.909  
23.844

24.779  25.714  26.649  27.584  28.519  29.454
zdef 7 levels
1000 925 850 700 500 300 200   === set Z 1 ou 2 ou 3 ... ou 7 (GRADS)
tdef 1 linear 6Z16may2009  6hr
vars 15
psnm  02,102,  0,  0 SEA LEVEL PRESSURE [hPa]
tsfc  0  187,  1,  0,  0 SURFACE TEMPERATURE [K]
prec  0   61,  1,  0,  0 TOTAL PRECIPITATION [Kg/m2/day]
cbnv  0   71,  3,  0,  0 CLOUD COVER [0-1]
role  0  114,  8,  0,  0 OUTGOING LONG WAVE AT TOP [W/m2]
mask  7  137,100 MASK [-/+]
uvel  7   33,100 ZONAL WIND (U) [m/s]
vvel  7   34,100 MERIDIONAL WIND (V) [m/s]
omeg  7   39,100 OMEGA [Pa/s]
fcor  7   35,100 STREAM FUNCTION [m2/s]
potv  7   36,100 VELOCITY POTENTIAL [m2/s]
zgeo  77,100 GEOPOTENTIAL HEIGHT [gpm]
temp  7   11,100 ABSOLUTE TEMPERATURE [K]
umrl  7   52,100 RELATIVE HUMIDITY [no Dim]
umes  7   51,100 SPECIFIC HUMIDITY [kg/kg]
endvars

Can someone help me to find out what is happening?

Tks a lot

Roberto Garcia


Frank Warmerdam wrote:

Roberto Garcia,MSc wrote:

Hi, tks for the answers.

According to my gdainfo output file attached, can anyone tell me what 
element could be a CLASSITEM?


Roberto,

In MapServer you always classify rasters based on the item [pixel] but
you do not need to explicitly state a classitem.

This is a relatively simple example of raster classification.

LAYER
  NAME grid1
  TYPE raster
  STATUS default
  DATA data/float.tif
  CLASS
NAME red
EXPRESSION ([pixel]  -3)
COLOR 255 0 0
  END
  CLASS
NAME green
EXPRESSION ([pixel] = -3 and [pixel]  3)
COLOR 0 255 0
  END
  CLASS
NAME blue
EXPRESSION ([pixel] = 3)
COLOR 0 0 255
  END
END

I can see my grib using Grads but what tool do I use to translate 
individual bands into GeoTIFF?


You can use gdal_translate to extract particular bands from a grib file
as a geotiff (or all the bands).  For instance, the following would
extract band 3 (total precipitation):

gdal_translate -b 3 T126L28.grb total_precip.tif

Converting the data on-the-fly to an image is the right way to work? 


Having MapServer access the grib file on the fly rather than depending
on pre-translated extracts will reduce data duplication and may help
avoid confusion about what is up to date.  However, it is possible that
there will be a performance penalty for extracting directly from GRIB.
If so, it may not be significant enough to matter to you.


What are the problems with Apache if I do not do so?


I don't see how either approach relates to Apache.  It is a
data management and possibly performance issue

Re: [mapserver-users] GRIB files

2009-05-13 Thread Roberto Garcia,MSc

Hi, tks for the answers.

According to my gdainfo output file attached, can anyone tell me what 
element could be a CLASSITEM?


I can see my grib using Grads but what tool do I use to translate 
individual bands into GeoTIFF?


Converting the data on-the-fly to an image is the right way to work? 
What are the problems with Apache if I do not do so? 


Tks

Frank Warmerdam wrote:

Jose Roberto M. Garcia, MSc wrote:

Hi all,

I'm trying to work with GRIB files in MapServer but it's not working. 
I know that when we work with images we need to classify the output 
using [pixel] as a CLASSITEM and make classes for each value but 
how is it for GRIB files? I have several levels an variables in it. 
How do I classify them? Does anybody explain me?


Roberto,

It would be helpful if you could make a GRIB file typical of what you 
want
to use available for download, so we could speak in the context of 
that file.


Grib files often include a number of variables.  These variables 
typically
get unrolled as bands.  If the variable is 3D (ie. has a time or 
elevation

dimension) then it may be unrolled as multiple bands from a GDAL point of
view.  One sample GRIB file I happen to have ends up with 40 bands, 
two of

which look like this when reported by gdalinfo:

Band 1 Block=360x1 Type=Float64, ColorInterp=Undefined
  Description = 5[m] DBSL (Depth below sea level)
  Metadata:
GRIB_UNIT=[K]
GRIB_COMMENT=Potential temperature [K]
GRIB_ELEMENT=POT
GRIB_SHORT_NAME=5-DBSL
GRIB_REF_TIME=  1130716800 sec UTC
GRIB_VALID_TIME=  1133308800 sec UTC
GRIB_FORECAST_SECONDS=2592000 sec
Band 2 Block=360x1 Type=Float64, ColorInterp=Undefined
  Description = 15[m] DBSL (Depth below sea level)
  Metadata:
GRIB_UNIT=[K]
GRIB_COMMENT=Potential temperature [K]
GRIB_ELEMENT=POT
GRIB_SHORT_NAME=15-DBSL
GRIB_REF_TIME=  1130716800 sec UTC
GRIB_VALID_TIME=  1133308800 sec UTC
GRIB_FORECAST_SECONDS=2592000 sec

As you can see these are temperature values for different depths
below sea level.

Once you know what band(s) you want to work with in MapServer you can
select them using the BANDS= processing directive.

  eg.
   PROCESSING BANDS=2

would select the 2nd band for display.

Then you can use classification as you describe.  In this case the
pixel values are presumably floating point temperatures measured in
degrees Kelvin though I'm not that familiar with the product.  Pay
particular attention to the online mapserver docs on classifying
non-8bit data as there are some caveats in this regard. I hope this
gets you going.

I would suggest:
  o inspecting your file(s) with gdalinfo, reviewing the metadata.
  o Translate individual bands into GeoTIFF files and visualize
and inspect the results with desktop software like QGIS or OpenEV.

Then move to mapserver once you are more comfortable you know how
the data is being represented through GDAL.

Best regards,
Driver: GRIB/GRIdded Binary (.grb)
Files: data\raster\T126L28.grb
Size is 161, 98
Coordinate System is:
GEOGCS[Coordinate System imported from GRIB file,
DATUM[unknown,
SPHEROID[Sphere,6367470,0]],
PRIMEM[Greenwich,0],
UNIT[degree,0.0174532925199433]]
Origin = (-150.000,-61.2460002)
Pixel Size = (0.938,-0.001836734693878)
Corner Coordinates:
Upper Left  (-150.000, -61.246) (150d 0'0.00W, 61d14'45.60S)
Lower Left  (-150.000, -61.426) (150d 0'0.00W, 61d25'33.60S)
Upper Right (   1.018, -61.246) (  1d 1'4.80E, 61d14'45.60S)
Lower Right (   1.018, -61.426) (  1d 1'4.80E, 61d25'33.60S)
Center  ( -74.491, -61.336) ( 74d29'27.60W, 61d20'9.60S)
Band 1 Block=161x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] MSL (Mean sea level)
  Metadata:
GRIB_UNIT=[Pa]
GRIB_COMMENT=Pressure reduced to MSL [Pa]
GRIB_ELEMENT=PRMSL
GRIB_SHORT_NAME=0-MSL
GRIB_REF_TIME=  1242453600 sec UTC
GRIB_VALID_TIME=  1242453600 sec UTC
GRIB_FORECAST_SECONDS=0 sec
Band 2 Block=161x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC (Ground or water surface)
  Metadata:
GRIB_UNIT=[-]
GRIB_COMMENT=undefined [-]
GRIB_ELEMENT=var187
GRIB_SHORT_NAME=0-SFC
GRIB_REF_TIME=  1242453600 sec UTC
GRIB_VALID_TIME=  1242453600 sec UTC
GRIB_FORECAST_SECONDS=0 sec
Band 3 Block=161x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] SFC (Ground or water surface)
  Metadata:
GRIB_UNIT=[kg/m^2]
GRIB_COMMENT=Total precipitation [kg/m^2]
GRIB_ELEMENT=APCP
GRIB_SHORT_NAME=0-SFC
GRIB_REF_TIME=  1242453600 sec UTC
GRIB_VALID_TIME=  1242453600 sec UTC
GRIB_FORECAST_SECONDS=0 sec
Band 4 Block=161x1 Type=Float64, ColorInterp=Undefined
  Description = 0[-] CTL (Cloud top level)
  Metadata:
GRIB_UNIT=[%]
GRIB_COMMENT=Total cloud cover [%]
GRIB_ELEMENT=TCDC
GRIB_SHORT_NAME=0-CTL
GRIB_REF_TIME=  1242453600 sec UTC
GRIB_VALID_TIME=  1242453600 sec UTC