Hello,

I have a TIFF or NetCD file that I want to read as a data set, do some 
work with it, and produce a new image (any format) as output. Can anyone 
suggest a simple library that could help me with that?

I'm running Ubuntu and I have access to the following libraries:

* libtiff-dev
* libtiff4
* GeoTIFF


I'm not familiar with any of these, and I don't even know if any of 
these have Fortran bindings. I also don't know exactly what a 
georeferenced TIFF file is, but I think that mine might be like that.


Let me explain a bit more about the problem I'm trying to solve:
I have this data set that describes the elevation or depth of every 
point on the Earth, which you can get from here:

http://ngdc.noaa.gov/mgg/global/global.html


My goal is to find out how much land and how much water Earth would have 
if it had more water or less water than it does today. This is part of a 
larger project where I try to estimate whether a typical Earth type 
planet would be a global ocean, a huge desert, or something in between 
like our Earth.

So basically I need to take this file with the elevation and depth and 
integrate the elevation data:


do i = min_elevation, max_elevation
     ... add the surface area with elevation between (i-1) and i  ...
end do


As a final bit, I'd like to draw a couple of sample maps of the Earth 
with different water contents. So image output is nice too, but less 
important.


As you can see, this task is in principle simple, except that I don't 
know what I'm doing. GIMP can't even open file file properly (16-bit 
data I assume) so it's hard to explore and feel my way through the problem.

In addition to recommending libraries, I would be happy to hear any 
advice that list members might have.

In case you are wondering: Yes, this is part of a school project. No, it 
does not violate any rules to ask for help. My project only requires a 
literature review, but I want to go the extra mile by adding a 
interesting programming component with an interesting idea.

Thank you for your help.

Cheers,
Daniel.

> Let me explain a bit more about the problem I'm trying to solve:
> I have this data set that describes the elevation or depth of every point 
> on the Earth, which you can get from here:

> My goal is to find out how much land and how much water Earth would have 
> if it had more water or less water than it does today. This is part of a 
> larger project where I try to estimate whether a typical Earth type planet 
> would be a global ocean, a huge desert, or something in between like our 
> Earth.

You might want to look for alternate sources of elevation data. IIRC the US 
Geological Survey has this type of data encoded in various ways to save 
space. So modern Fortran really helps with "stream" I/O and bit-field 
manipulation.

-- e

On 12/16/2011 03:45 PM, e p chandler wrote:
> You might want to look for alternate sources of elevation data. IIRC the US
> Geological Survey has this type of data encoded in various ways to save
> space. So modern Fortran really helps with "stream" I/O and bit-field
> manipulation.

Is there a reason why a different format would work better? Is there a 
format that is particularly easy to import from Fortran? The NOAA 
website has netCDF, GRD98, "Georeferenced TIFF" and "xyz" (whatever that 
means).

Daniel.

On 12/16/2011 05:12 PM, Daniel Carrera wrote:
> On 12/16/2011 03:45 PM, e p chandler wrote:
>> You might want to look for alternate sources of elevation data. IIRC
>> the US
>> Geological Survey has this type of data encoded in various ways to save
>> space. So modern Fortran really helps with "stream" I/O and bit-field
>> manipulation.
>
> Is there a reason why a different format would work better? Is there a
> format that is particularly easy to import from Fortran? The NOAA
> website has netCDF, GRD98, "Georeferenced TIFF" and "xyz" (whatever that
> means).

I think I'm going to try the NetCDF format. Ubuntu has some packages to 
manage NetCDF files including a library for Fortran 90 and one for Python.

Cheers,
Daniel.

The ETOPO1 data are also available in binary form: just open the tap
and let the bytes flow! No need for a library! I usually generate .bmp-
files with my own subroutines, these files are quite easy in format.
For more elaborate plots I generate Postscript files and convert them
to jpg with a system-call to utility "gs" afterwards. Make sure you
have the "bedrock" flavour of the elevation data, otherwise you'll
remain stuck on top of today's ice sheets when you lower the water
level! In case you want the .bmp- and Postscript-routines, please let
me know!

Arjan

On 12/16/2011 06:56 PM, Arjan wrote:
> The ETOPO1 data are also available in binary form: just open the tap
> and let the bytes flow! No need for a library!

Sounds good, but how do I read that? This is probably a dumb question 
but this is the first time I do anything remotely like this.

> I usually generate .bmp-
> files with my own subroutines, these files are quite easy in format.
> For more elaborate plots I generate Postscript files and convert them
> to jpg with a system-call to utility "gs" afterwards.

Ok, thanks. Both options sound reasonable.


> Make sure you have the "bedrock" flavour of the elevation data, otherwise
> you'll remain stuck on top of today's ice sheets when you lower the water
> level!

Yes, very important.

> In case you want the .bmp- and Postscript-routines, please let
> me know!

Yes, I'd like to have those. Thanks!

Daniel.

> Sounds good, but how do I read that? This is probably a dumb question
> but this is the first time I do anything remotely like this.

I guess that if you define a 2D array of REAL(4) that will fit this:

ncols         21601
nrows         10801
xllcenter     -180.00
yllcenter     -90.00
cellsize      0.01666666667
NODATA_value  -99999
byteorder     LSBFIRST
NUMBERTYPE    4_BYTE_FLOAT
MIN_VALUE     -10898.0
MAX_VALUE     8271.0

then a single read-statement will do the work. Check the order of the
rows and columns! Latitude seems to run from North Pole to Antarctica!


> > Make sure you have the "bedrock" flavour of the elevation data, otherwise
> > you'll remain stuck on top of today's ice sheets when you lower the water
> > level!
>
> Yes, very important.

As a dutchman, I have to ask your attention for the dry parts that can
be found below sea-level! Once you'll have successfully read the
elevation data, may I ask you to do a test for the Netherlands?

> > In case you want the .bmp- and Postscript-routines, please let
> > me know!
>
> Yes, I'd like to have those. Thanks!

I have to polish them up a bit to make them usable for other people
than me ;-) . If you start to prepare/make your calculations, I'll
compile a set of sources plus an example.

Regards,


Arjan

Hi Arjan,

Thanks for the help. I've looked at the ".xyz" file, which is a big 
ASCII file and I can see how the data is organized:

-179.992	89.9917	-4229
-179.975	89.9917	-4229
-179.958	89.9917	-4229


The columns are clearly longitude (-180 to 180 deg), latitude (90 to 90 
deg) and I've figured that the height is in meters relative to sea 
level. So I just have one question right now:


> I guess that if you define a 2D array of REAL(4) that will fit this:
>
> ncols         21601
> nrows         10801
> xllcenter     -180.00
> yllcenter     -90.00
> cellsize      0.01666666667


What is the unit of "cellsize"? Do you think the cells are all the same 
area or the same angular size? I'm wondering if I have to compute the 
area of each cell.


> As a dutchman, I have to ask your attention for the dry parts that can
> be found below sea-level! Once you'll have successfully read the
> elevation data, may I ask you to do a test for the Netherlands?

A slightly unrelated question:  Can you explain to me why the 
Netherlands is not a huge lake? Since it is below sea level, I would 
have imagined that the moment it rains the water would get stuck there 
and eventually it'd form a huge lake.

Cheers,
Daniel.

"Daniel Carrera"  wrote in message

>> As a dutchman, I have to ask your attention for the dry parts that can
>> be found below sea-level! Once you'll have successfully read the
>> elevation data, may I ask you to do a test for the Netherlands?

>A slightly unrelated question:  Can you explain to me why the 
>Netherlands is not a huge lake? Since it is below sea level, I would 
>have imagined that the moment it rains the water would get stuck there 
>and eventually it'd form a huge lake.

That's what windmills are for! [smile]

Ok, I think I can answer my own question:

The cells are 0.017 degrees in width and height, so yes, I have to 
compute the surface area of each cell. That means that the are of each 
cell ranges between 0.0142 m^2 at the equator to 0.00000206 m^2 at the 
poles.

Wow, that's far more precise than anything I need, but oh well, this is 
the data I have. I'd rather the data be to precise than not precise enough.

Cheers,
Daniel.


On 12/18/2011 10:05 PM, Daniel Carrera wrote:
> Hi Arjan,
>
> Thanks for the help. I've looked at the ".xyz" file, which is a big
> ASCII file and I can see how the data is organized:
>
> -179.992 89.9917 -4229
> -179.975 89.9917 -4229
> -179.958 89.9917 -4229
>
>
> The columns are clearly longitude (-180 to 180 deg), latitude (90 to 90
> deg) and I've figured that the height is in meters relative to sea
> level. So I just have one question right now:
>
>
>> I guess that if you define a 2D array of REAL(4) that will fit this:
>>
>> ncols 21601
>> nrows 10801
>> xllcenter -180.00
>> yllcenter -90.00
>> cellsize 0.01666666667
>
>
> What is the unit of "cellsize"? Do you think the cells are all the same
> area or the same angular size? I'm wondering if I have to compute the
> area of each cell.
>
>
>> As a dutchman, I have to ask your attention for the dry parts that can
>> be found below sea-level! Once you'll have successfully read the
>> elevation data, may I ask you to do a test for the Netherlands?
>
> A slightly unrelated question: Can you explain to me why the Netherlands
> is not a huge lake? Since it is below sea level, I would have imagined
> that the moment it rains the water would get stuck there and eventually
> it'd form a huge lake.
>
> Cheers,
> Daniel.

Ok, I can't multiply.

After realizing how absurd the values I posted were, I took a second 
look and found my error. I now think that the area of each cell ranges 
from 3.57 km^2 at the equator to 517 m^2 at the poles. That sounds a lot 
more reasonable now.

Cheers,
Daniel.


On 12/18/2011 10:14 PM, Daniel Carrera wrote:
> Ok, I think I can answer my own question:
>
> The cells are 0.017 degrees in width and height, so yes, I have to
> compute the surface area of each cell. That means that the are of each
> cell ranges between 0.0142 m^2 at the equator to 0.00000206 m^2 at the
> poles.
>
> Wow, that's far more precise than anything I need, but oh well, this is
> the data I have. I'd rather the data be to precise than not precise enough.
>
> Cheers,
> Daniel.
>
>
> On 12/18/2011 10:05 PM, Daniel Carrera wrote:
>> Hi Arjan,
>>
>> Thanks for the help. I've looked at the ".xyz" file, which is a big
>> ASCII file and I can see how the data is organized:
>>
>> -179.992 89.9917 -4229
>> -179.975 89.9917 -4229
>> -179.958 89.9917 -4229
>>
>>
>> The columns are clearly longitude (-180 to 180 deg), latitude (90 to 90
>> deg) and I've figured that the height is in meters relative to sea
>> level. So I just have one question right now:
>>
>>
>>> I guess that if you define a 2D array of REAL(4) that will fit this:
>>>
>>> ncols 21601
>>> nrows 10801
>>> xllcenter -180.00
>>> yllcenter -90.00
>>> cellsize 0.01666666667
>>
>>
>> What is the unit of "cellsize"? Do you think the cells are all the same
>> area or the same angular size? I'm wondering if I have to compute the
>> area of each cell.
>>
>>
>>> As a dutchman, I have to ask your attention for the dry parts that can
>>> be found below sea-level! Once you'll have successfully read the
>>> elevation data, may I ask you to do a test for the Netherlands?
>>
>> A slightly unrelated question: Can you explain to me why the Netherlands
>> is not a huge lake? Since it is below sea level, I would have imagined
>> that the moment it rains the water would get stuck there and eventually
>> it'd form a huge lake.
>>
>> Cheers,
>> Daniel.
>

Please also download the ETOPO5 dataset: quite a bit smaller, but it
will allow you to make fast test-runs!

Concerning the area of your tiles: R^2 x COS(latitude) is your friend,
R being the earth's radius.
This assumes that you measure your angles in radians and not in
degree.

E.P Chandler is right about why we can live at 10 meters below sea.
With sea-level rise as a possible scenario for the near future, I
would not be surprised to see some major changes come...

Arjan

Arjan <arjan.van.dijk@rivm.nl> wrote:

(snip)
> Concerning the area of your tiles: R^2 x COS(latitude) is your friend,
> R being the earth's radius.
> This assumes that you measure your angles in radians and not in
> degree.

There is a correction for the oblate spheroidal shape, though
maybe not needed for this problem.  

Otherwise, you might want to consider them as trapazoids, with 
cos(lattitude at the bottom) and cos(lattitude at the top).
A little closer than cos(center lattitude)

-- glen

> There is a correction for the oblate spheroidal shape, though
> maybe not needed for this problem.
>
> Otherwise, you might want to consider them as trapazoids, with
> cos(lattitude at the bottom) and cos(lattitude at the top).
> A little closer than cos(center lattitude)


I get the impression that a rough question like "What if the Earth
would have received twice as much water?" does not require an answer
with more than 1 or 2 significant digits. Therefore, there is no need
for a HUGE dataset like ETOPO1 either! ETOPO5 (resolution already
better than 10 km!) will do well enough.

A.

On 12/20/2011 12:08 PM, Arjan wrote:
>> Otherwise, you might want to consider them as trapazoids, with
>> cos(lattitude at the bottom) and cos(lattitude at the top).
>> A little closer than cos(center lattitude)
>
> I get the impression that a rough question like "What if the Earth
> would have received twice as much water?" does not require an answer
> with more than 1 or 2 significant digits. Therefore, there is no need
> for a HUGE dataset like ETOPO1 either! ETOPO5 (resolution already
> better than 10 km!) will do well enough.


Indeed. There is no need for high precision here.

There is some research trying to estimate how much water a terrestrial 
planet might receive. I think for Earth the ballpark answer is between 
0.5 and 3 times the water that it has today, but it's all very rough.

For my project I only need to read the relevant papers and summarize 
them, but I want to go a bit further and explain what it would mean for 
Earth to have 0.5 or 3 times the water it has today.

Cheers,
Daniel.

> them, but I want to go a bit further and explain what it would mean for
> Earth to have 0.5 or 3 times the water it has today.


Daniel, from the sample code (based on ETOPO5) that I sent you, you
can see that:

0.5 x the present amount of water would require a sealevel 2090 m
below today's level, covering 58.8% of the Earth's surface (71.1% now)

3 x  the present amount of water would require a sealevel 5550 m above
today's level, covering 98.8% of the Earth's surface (71.1% now). The
fraction of 98.8% may actually be higher, since somehow my code does
not reach 100% for nearly 9 km of sealevel rise...

The sample code can also be used to construct maps of the Earth for
these 2 conditions! Please show me the resulting maps once you'll
have'm! I'm interested to see what they look like!


Arjan

ps. In reality, the thermo-haline cycle has lead to a distribution of
salt and heat such that re-distribution would lead to a substantial
sea-level rise without the need to add water to our planet!

Arjan