Pyresample to geotiff

Rutger

unread,

Mar 5, 2012, 8:21:59 AM3/5/12

to pyt...@googlegroups.com

Hello,

I successfully managed to use Pyresample to convert swath data stored in a .csv file (from a HDF-table) to a nice plot using matplotlib (+basemap), where i can save it as a PNG file.
However, besides a nice visual plot i would also like to save the grid in some format which supports geolocation, like Geotiff. Is this possible from Pyresample? I couldnt find anything on saving grids to a file. Or should i move my 'result' from Pyresample to something in mpop?

My script currently looks like:

import numpy as np
import pyresample as pr
import matplotlib.pyplot as plt

# read data
csvfile = 'D:\\swath.csv'
data = np.genfromtxt(csvfile, delimiter=',', dtype=None, names=True)

lons = data['Longitude']
lats = data['Latitude']
data = data['Soil_moisture']

# read area definition (i modified the default ease_nh so it is 1700*1700)
area_def = pr.utils.load_area('D:\\areas.cfg', 'ease_nh_hires')

# reduce data to area
grid_lons, grid_lats = area_def.get_lonlats()
red_lons, red_lats, red_data = pr.data_reduce.swath_from_lonlat_grid(grid_lons, grid_lats, lons, lats, data, radius_of_influence=12500)

# project
swath_def = pr.geometry.SwathDefinition(lons=red_lons, lats=red_lats)
result = pr.kd_tree.resample_gauss(swath_def, red_data, area_def, radius_of_influence=25000, sigmas=12500, fill_value=None)

# plot quicklook
fig1 = plt.figure(1, figsize=(3200 / 200, 2400 / 200))
bmap = pr.plot.area_def2basemap(area_def)
bmng = bmap.bluemarble()
col = bmap.imshow(result, origin='upper', vmin=0, vmax=100)
plt.title('Pyresample test')
plt.colorbar()
# plt.show()
plt.savefig('test.png', dpi=200, bbox_inches='tight')

All i can think of is constructing the geolocated file in something like GDAL and force it on my numpy result from "pr.kd_tree.resample_gauss", but would be a rather ugly workaround.

Any help would be appreciated, thanks,
Rutger

Martin Raspaud

unread,

Mar 5, 2012, 9:01:57 AM3/5/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 05/03/12 14:21, Rutger wrote:
> Hello,

Hi Rutger,

> I successfully managed to use Pyresample to convert swath data
> stored in a .csv file (from a HDF-table) to a nice plot using
> matplotlib (+basemap), where i can save it as a PNG file. However,
> besides a nice visual plot i would also like to save the grid in
> some format which supports geolocation, like Geotiff. Is this
> possible from Pyresample? I couldnt find anything on saving grids
> to a file. Or should i move my 'result' from Pyresample to
> something in mpop?

Yes, I think you need mpop to save to geotiff. The quickest way would
be to create a GeoImage object, inject the data in it, and save it as
geotiff. See the geo_image interface here:
http://mraspaud.github.com/mpop/image.html#module-mpop.imageo.geo_image

Otherwise, the nice way would be to add a reader for your data to
mpop. Some indications here:
http://pytroll.org/quickstart_custom.html

Keep us updated on how it goes !

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVMdVAAoJEBdvyODiyJI4dQgH+wWXhBA0sf+c75Aq4KW19Gig
j1nk3tBifMCX1PA8mKjDzdTcDIjD54hQgSvrC7hTXila7LXu99KLyKfWmg85GINn
6pjbR1vcTugY1Xc0ieXytZ0e+CkMm82JRZhZlBDkt3MnVM1EIzfOMwlfKQ4m0CZG
J7giIRlOlZUNMXhJ9ZtmHylzMnLNx1cOTgg62Hw5Y5U7trB+BWPGvCJIxxGAupRd
L/vQJLzU2DGAaeOjOE0ghXnXyEBgKPvc6UW9m4n4wBv9Xqju2bxZzV+6hFfpH4rv
hlG7yiUw6+23bP5H9X5oHWYNFe9i0vxYlHJ4NxlgfYY/qqzw+33jTYhY59ZnP5M=
=Ki1H
-----END PGP SIGNATURE-----

martin_raspaud.vcf

Rutger

unread,

Mar 5, 2012, 10:22:45 AM3/5/12

to pyt...@googlegroups.com

Hi Martin,

Thank you for your response. I am already stuck at the installation of mpop. Here is what i have done so far:

I'm on Windows using Python 2.7.

I have downloaded the mpop files using the "download as .tar.gz" at bithub. After extracting i ran the setup.py file which displays this error:
"ValueError: Cannot find the version number!"

After looking in the version.py this error makes sense. (I guess) Since i'm working in an unpacked 'tarball' instead of a 'git working copy' the 'call_git_describe' function fails and returns None. It should then fallback on the version mentioned in RELEASE-VERSION, but this file didn't come with the package, so it fails as well. To overcome this i tried manually creating this file containing only 1 line with '0.10.0', running 'python setup.py --version' then picks up this version number. After creating this file running "python setup.py install" gets past this version checking and starts installing. After a while it displays this error:

installing package data to build\bdist.win32\egg
running install_data
creating build\bdist.win32\egg\etc
copying etc\geo_image.cfg -> build\bdist.win32\egg\etc
copying etc\world_map.ascii -> build\bdist.win32\egg\etc
Traceback (most recent call last):
File "setup.py", line 68, in <module>
    'mipp']
File "C:\Python27\lib\distutils\core.py", line 152, in setup
    dist.run_commands()
File "C:\Python27\lib\distutils\dist.py", line 953, in run_commands
    self.run_command(cmd)
File "C:\Python27\lib\distutils\dist.py", line 972, in run_command
    cmd_obj.run()
File "C:\Python27\Lib\site-packages\setuptools\command\install.py", line 76, in run
    self.do_egg_install()
File "C:\Python27\Lib\site-packages\setuptools\command\install.py", line 96, in do_egg_install
    self.run_command('bdist_egg')
File "C:\Python27\lib\distutils\cmd.py", line 326, in run_command
    self.distribution.run_command(command)
File "C:\Python27\lib\distutils\dist.py", line 972, in run_command
    cmd_obj.run()
File "C:\Python27\Lib\site-packages\setuptools\command\bdist_egg.py", line 195, in run
    self.do_install_data()
File "C:\Python27\Lib\site-packages\setuptools\command\bdist_egg.py", line 145, in do_install_data
    self.call_command('install_data', force=0, root=None)
File "C:\Python27\Lib\site-packages\setuptools\command\bdist_egg.py", line 161, in call_command
    self.run_command(cmdname)
File "C:\Python27\lib\distutils\cmd.py", line 326, in run_command
    self.distribution.run_command(command)
File "C:\Python27\lib\distutils\dist.py", line 972, in run_command
    cmd_obj.run()
File "C:\Python27\lib\distutils\command\install_data.py", line 58, in run
    dir = convert_path(f[0])
File "C:\Python27\lib\distutils\util.py", line 201, in convert_path
    raise ValueError, "path '%s' cannot end with '/'" % pathname
ValueError: path 'share/doc/mpop/' cannot end with '/'

Im not sure what to do with this error. In the setup.py i see this ('share/doc/'+NAME+'/', at line 57, maybe the path created isnt correct? Do you have any idea.

fyi: running the setup.py from Pyresample worked 'out of the box'.

Rutger

Op maandag 5 maart 2012 15:01:57 UTC+1 schreef Martin Raspaud het volgende:

Op maandag 5 maart 2012 15:01:57 UTC+1 schreef Martin Raspaud het volgende:

Martin Raspaud

unread,

Mar 5, 2012, 10:48:37 AM3/5/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 05/03/12 16:22, Rutger wrote:
> Hi Martin,

Hi Rutger,

> Thank you for your response. I am already stuck at the
> installation of mpop. Here is what i have done so far:
>
> I'm on Windows using Python 2.7.

Ok, I have to say mpop is totally untested on windows, so I'm glad to
get your mail :o)

> I have downloaded the mpop files using the "download as .tar.gz"
> at bithub. After extracting i ran the setup.py file which displays
> this error: "ValueError: Cannot find the version number!"
>
> After looking in the version.py this error makes sense. (I guess)
> Since i'm working in an unpacked 'tarball' instead of a 'git
> working copy' the 'call_git_describe' function fails and returns
> None. It should then fallback on the version mentioned in
> RELEASE-VERSION, but this file didn't come with the package, so it
> fails as well. To overcome this i tried manually creating this
> file containing only 1 line with '0.10.0', running 'python
> setup.py --version' then picks up this version number. After
> creating this file running "python setup.py install" gets past this
> version checking and starts installing.

Ok, the version number thing is indeed counting on a git cloning,
which is apparently not optimal, I'll change this.

> After a while it displays this error:
>
> installing package data to build\bdist.win32\egg

[...]

>> File "C:\Python27\lib\distutils\util.py", line 201, in
>> convert_path raise ValueError, "path '%s' cannot end with '/'" %
>> pathname ValueError: path 'share/doc/mpop/' cannot end with '/'
>>
>
> Im not sure what to do with this error. In the setup.py i see this
> ('share/doc/'+NAME+'/', at line 57, maybe the path created isnt
> correct? Do you have any idea.

Yes, I put "/" which is not working in windows of course.
I should replace all these by os.path.join.
I attach a new setup.py, can you try using it ?

> fyi: running the setup.py from Pyresample worked 'out of the box'.

Yes, I'm doing some exotic stuff in the mpop setup for the
documentation and config files there...

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVOBVAAoJEBdvyODiyJI4SCAIAJ4oVaylsNuPox0fT7GIOs0J
WYdU7Ne266qZUREhe5lUw0HKZ82/YG4jRYHz4l7n2uSjFpyW1y0dbm9QkPLRFwcI
8JI8z94XBPZZvDkjAiDhzR4ZuxqGumgGof+u+Izb0yg9kho5jPYq94AqMKk2egb+
EQX28e3OYkITnJb7H2abbYcLNE9Dzp8+LykXHx3fr30aCXhvDkmtmysHbm33b39o
+jhW11YhKI+Y3vY0eNjRp0+dK6m+lNQDwUcg2eeRVlaV8RP3M6WFNl2vSxTiKzjB
MYeNtgiez/iVBnhCBNyZuXXI4VvUYlpw8vXMK0JySx+UA1L3PIb12morgG2shjI=
=dMm0
-----END PGP SIGNATURE-----

setup.py

setup.py.sig

martin_raspaud.vcf

Rutger

unread,

Mar 5, 2012, 11:49:30 AM3/5/12

to pyt...@googlegroups.com

Hey,

Thank you for your fast response. The new setup.py did the trick, installation went smooth. Importing mpop on a 'blank' python command line gave no errors, so that seems alright.

As you can see in my initially posted script i import matplotlib (import matplotlib as plt). When i added 'import mpop' to the script i got this error:

"AccessInit: hash collision: 3 for both 1 and 1"

Apparently, some of the imports are conflicting, i found that commenting out the matplotlib line solved this error (although i would like to use it eventually).

I then inserted the following two lines right after my 'result = pr.kd_tree.resample_gauss(....':

img = GeoImage(result, 'ease_nh_hires', time_slot=None, mode='L', crange=None, fill_value=None, palette=None)
img.geotiff_save('test.tif', compression=6, tags=None, gdal_options=None, blocksize=0, geotransform=None, spatialref=None)

Resulting in this error mesage:

[WARNING: 2012-03-05 17:28:31 : mpop] Couldn't find the mpop.cfg file. Do you have one ? is it in $PPP_CONFIG_DIR ?

Traceback (most recent call last):

File "test.py", line 31, in <module>
    img.geotiff_save('test.tif', compression=6, tags=None, gdal_options=None, blocksize=0, geotransform=None, spatialref=None)
File "C:\Python27\lib\site-packages\mpop-v0.10.0-py2.7.egg\mpop\imageo\geo_image.py", line 234, in geotiff_save
    area = get_area_def(self.area_id)
File "C:\Python27\lib\site-packages\mpop-v0.10.0-py2.7.egg\mpop\projector.py", line 62, in get_area_def
    return utils.parse_area_file(AREA_FILE, area_name)[0]
File "C:\Python27\lib\site-packages\pyresample-0.7.10-py2.7.egg\pyresample\utils.py", line 75, in parse_area_file
    area_file = open(area_file_name, 'r')
IOError: [Errno 22] invalid mode ('r') or filename: ''

I cant find the mpop.cfg myself either. I dont have a PPP_CONFIG_DIR as an environment variable, so it should look in the etc\ folder of the installation if im correct, which it isnt. Im not sure if the parameters supplied to GeoImage are correct, especially the area_id, but it doesnt throw an error at this line. And the error above hints at something else.

Anyway, i wont be able to do any further tests today, thank you for your assistance so far.

Rutger

Martin Raspaud

unread,

Mar 6, 2012, 1:01:31 AM3/6/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 05/03/12 17:49, Rutger wrote:
> Hey,

Hi Rutger,

> Thank you for your fast response. The new setup.py did the trick,
> installation went smooth. Importing mpop on a 'blank' python
> command line gave no errors, so that seems alright.
>
> As you can see in my initially posted script i import matplotlib
> (import matplotlib as plt). When i added 'import mpop' to the
> script i got this error:
>
>> "AccessInit: hash collision: 3 for both 1 and 1"
>>
>
> Apparently, some of the imports are conflicting, i found that
> commenting out the matplotlib line solved this error (although i
> would like to use it eventually).

Strange, I can't reproduce this error...

> I then inserted the following two lines right after my 'result =
> pr.kd_tree.resample_gauss(....':
>
>> img = GeoImage(result, 'ease_nh_hires', time_slot=None,
>> mode='L', crange=None, fill_value=None, palette=None)
>> img.geotiff_save('test.tif', compression=6, tags=None,
>> gdal_options=None, blocksize=0, geotransform=None,
>> spatialref=None)
>>
>
> Resulting in this error mesage:
>
>> [WARNING: 2012-03-05 17:28:31 : mpop] Couldn't find the mpop.cfg
>> file. Do you have one ? is it in $PPP_CONFIG_DIR ?

[...]

>> File
>> "C:\Python27\lib\site-packages\pyresample-0.7.10-py2.7.egg\pyresample\utils.py",
>>
>>
line 75, in parse_area_file area_file = open(area_file_name,
>> 'r') IOError: [Errno 22] invalid mode ('r') or filename: ''
>>
>
>
> I cant find the mpop.cfg myself either. I dont have a
> PPP_CONFIG_DIR as an environment variable, so it should look in
> the etc\ folder of the installation if im correct, which it isnt.
> Im not sure if the parameters supplied to GeoImage are correct,
> especially the area_id, but it doesnt throw an error at this line.
> And the error above hints at something else.

Ok, the mpop.cfg file is not provided as is, but there is a
mpop.cfg.template you can copy in the etc directory.
As you will see, it says in which file your areas are defined, by
default areas.def (which is only provided as a template too). The area
you specify in the GeoImage constructor has to be defined there. This
is if you want you area to be stored somewhere.

However, a way to avoid all these files is if you provide the area
(not only the string of the area name):

>>> from pyresample import geometry area_id = 'ease_sh' name =
>>> 'Antarctic EASE grid' proj_id = 'ease_sh' proj4_args =
>>> 'proj=laea, lat_0=-90, lon_0=0, a=6371228.0, units=m' x_size =
>>> 425 y_size = 425 area_extent =
>>> (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> proj_dict = {'a': '6371228.0', 'units': 'm', 'lon_0': '0',
... 'proj': 'laea', 'lat_0': '-90'}
>>> area_def = geometry.AreaDefinition(area_id, name, proj_id,
proj_dict, x_size,
... y_size, area_extent)
>>> print area_def
Area ID: ease_sh
Name: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'a': '6371228.0', 'units': 'm', 'lon_0': '0', 'proj':
'laea', 'lat_0': '-90'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

And from there, you can just run:
>>> img = GeoImage(result, area_def, time_slot=None, mode='L')
>>> img.geotiff_save('test.tif', compression=6)

Hope that helps !

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVag7AAoJEBdvyODiyJI4uT4IAIn6pP8t055ige+rmNEf3qbL
botP4TCT9w8eZZI4hFMrYdvsTWjnCrDmfOZ9nlT2ZS1U+SqXmLvCV9RnIIak/Hmp
7WRSY9nY/GPpoY7vUqRsowW5byYcQhEE/bgJ0mOwHqqrMT48dXdi19//XRAYgAI9
fuIzU+6GcwsvI/CbyLmWLQoyNwN3dkw2GAgrmk3w+kL25Bf8TeHeNA0IK4fFKYP0
tJiNLk2ngeoGoyqwslCYSfRU5L7vcMQKxtdjKh/yByLGlYFolCH1ZscEJjRs/6Mg
7m8untbMRUE3W4SxBSKi56B7DsthFqlKgM1bQw7iyy5wPdPxhOF4DRk7llJexzk=
=zNvf
-----END PGP SIGNATURE-----

martin_raspaud.vcf

Martin Raspaud

unread,

Mar 6, 2012, 1:01:31 AM3/6/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 05/03/12 17:49, Rutger wrote:

Hi Rutger,

Strange, I can't reproduce this error...

[...]

Ok, the mpop.cfg file is not provided as is, but there is a
mpop.cfg.template you can copy in the etc directory.
As you will see, it says in which file your areas are defined, by
default areas.def (which is only provided as a template too). The area
you specify in the GeoImage constructor has to be defined there. This
is if you want you area to be stored somewhere.

However, a way to avoid all these files is if you provide the area
(not only the string of the area name):

... 'proj': 'laea', 'lat_0': '-90'}

proj_dict, x_size,
... y_size, area_extent)

Area ID: ease_sh

Name: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'a': '6371228.0', 'units': 'm', 'lon_0': '0', 'proj':
'laea', 'lat_0': '-90'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

And from there, you can just run:

Hope that helps !

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVagzAAoJEBdvyODiyJI4/2oIAMDsXJgAS8jA9ejaf0hE+XsQ
gtRmSqLP0mvadusS3cD6mHY/bYWKmiLYFz5UDMrr5/Y5vHVtdxa67LoeUqwAkQJm
MyKjhG0TV0HBSd920Rtfev+sF9A0JRhjGdTQX+DwPU0tbUqkzCsCFOpwr6mJPNL/
h36bRTc9Qq5Y5xGgCyby4glREMiakgEMOm3Kf03fVgv+eCMMoVQwz/r5CYbRLsVj
xboyE8PNjwEaBP6zkaoeri9cH/VOz+eOI5dw35V52RuUOIHxxr0zt+cpw2YD9XIC
/QA0cZwEkHgyOQedAbYxg3B9yoGSFj1ZLVMi3recGwwjSZ2WWF2VQuEUznS0vmk=
=JgJn
-----END PGP SIGNATURE-----

martin_raspaud.vcf

Martin Raspaud

unread,

Mar 6, 2012, 4:04:25 AM3/6/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi again,

I see now that googlegroups destroyed the layout of the python code...

I actually copied the area definition from the pyresample documentation:
http://pyresample.googlecode.com/git/docs/build/html/geo_def.html

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVdMZAAoJEBdvyODiyJI4dbkIANpMs2em8EF2uIAcrpsppqQM
YU4eEHykbYYXOBpjS3FRIUANLHNpuWfFzDdgSckSD1tvSEMjH5QU/XrrK4BFO9x+
utygzLN1fzNiAsN1Uc0rgeLaBvPnrLeKE2AAGuX5klhV9j4bRWFDsHlc8kBJ1XXd
Wugj+fzb64iaM6S/aoh9z8GSNTgcmI8veXNwbK1qML1n3/PQciyiPHzvVW7bWSy2
vcs23Z6nsfsBvDQKrBA9BmkjDJC8rdy4aacbkjDL+2W4vg+l5pqZTFsgZzpC7VT6
LVRhPHEcTbtnrXOX+2L+EQRZpFLMMnpz+WSOWWbRBMsfb6+xeU9ioc3by9R9AdM=
=4D8j
-----END PGP SIGNATURE-----

martin_raspaud.vcf

Rutger

unread,

Mar 6, 2012, 4:36:57 AM3/6/12

to pyt...@googlegroups.com

Hey,

I will see if the matplotlib conlict has something to do with my current Python environment. The worst case would be splitting up script in a quicklook and a geotiff part, which would work for me.

Defining a proper 'area definition' resolved the previous issue about not finding mpop.cfg. Only after your post i realized what i was doing wrong, i now reused my area definition which was already defined for the initial swath resampling. You want it to be equal so redefining it would only make thinks worse.

Using "time_slot=None" results in "AttributeError: 'NoneType' object has no attribute 'strftime'", so i just provided a datetime object with dummy values, this works.

The output is now a proper Geotiff although the projection isnt recognized by any other software im using. I tried gdalinfo, uDig and Erdas Imagine. This might be caused by the definition of the NH EASE-grid projection, i will try some other projections later. I only used the NH EASE-grid because it was used in the Pyresample documentation. For me its not really a problem because software like GDAL easily allows overriding the projection attached to a file.

One last thing, and im afraid this is by design. My input data has integer values between 0 and 100, to get exactly these values in my geotiff i scaled them using crange=[0,255]. But this trick only works if your input data consist of integers lower then 255. Is there a way to save a Geotiff with for example floating values, or 16bit 'digital numbers' etc? And i would also like to disable the alpha band, which is now added by default i think?

I feel like im abusing the GeoImage object for things its not designed for. :) Maybe after resamling my raw swath to a numpy grid swithing to GDAL will be easier after all, what do you think?

Again, thanks for your support.

Rutger

Op dinsdag 6 maart 2012 07:01:31 UTC+1 schreef Martin Raspaud het volgende:

Op dinsdag 6 maart 2012 07:01:31 UTC+1 schreef Martin Raspaud het volgende:

Martin Raspaud

unread,

Mar 6, 2012, 5:20:42 AM3/6/12

to pyt...@googlegroups.com

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 06/03/12 10:36, Rutger wrote:
> Hey,

Hi Rutger,

[...]

> Using "time_slot=None" results in "AttributeError: 'NoneType'
> object has no attribute 'strftime'", so i just provided a datetime
> object with dummy values, this works.

Oh, that could be fixed I guess :)

> The output is now a proper Geotiff although the projection isnt
> recognized by any other software im using. I tried gdalinfo, uDig
> and Erdas Imagine. This might be caused by the definition of the NH
> EASE-grid projection, i will try some other projections later. I
> only used the NH EASE-grid because it was used in the Pyresample
> documentation. For me its not really a problem because software
> like GDAL easily allows overriding the projection attached to a
> file.

We use these geotiff against a mapserver (and other visualisation
software), and it's working so far. But we use a polar stereographic
projection...

> One last thing, and im afraid this is by design. My input data has
> integer values between 0 and 100, to get exactly these values in my
> geotiff i scaled them using crange=[0,255]. But this trick only
> works if your input data consist of integers lower then 255. Is
> there a way to save a Geotiff with for example floating values, or
> 16bit 'digital numbers' etc?

It might not be so hard to add an option to the geotiff save function
to support floating values... let's see. Try with the attached file,
and img.save("myimage.tif", floating_point=True) for example.

> And i would also like to disable the alpha band, which is now
> added by default i think?

The alpha band can be disabled by providing a fill_value that is not
None. Try putting fill_value to (0.0, ) for a grayscale image.

> I feel like im abusing the GeoImage object for things its not
> designed for. :) Maybe after resamling my raw swath to a numpy grid
> swithing to GDAL will be easier after all, what do you think?

No, not abusing... Floating point is supported by geotiff, so there's
no reason it shouldn't be supported by geo_image :)

> Again, thanks for your support.

No problem, thanks for you feedback, it's highly valuable to us !

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJPVeT5AAoJEBdvyODiyJI4bfwIANm505fQj5VpscBdGRLHDCDs
senaWsguSfAYazjQVXeZsfxvFW0fCRubum80xBt+/RlsJdnQGBup9ZHH5kac7OAU
b5rZSmfR8JohfRpSDdJ8ViGMfKEGkn9XUTBhpnDmRWk6+iwiqYEo/aNzwDxFkX7+
UlfzYBsBaPNyuNjjtwlSuZZanCEeNNmtpRQxMwgZJei9DaFIGPT1hFmF/mnjKajK
qWMS90v44F05X04apGxerq67QHgnNGLQspyEj36MwJpsZBt1rWasbQRcUclbEaYW
vmfaHDeQE3NbzRXO5hEEPw52xMmts9hW3qIZRUG7gVgaD0AHW/lgXzoZESiTJMI=
=J3dr
-----END PGP SIGNATURE-----

geo_image.py

geo_image.py.sig

martin_raspaud.vcf

Rutger

unread,

Mar 6, 2012, 7:07:46 AM3/6/12

to pyt...@googlegroups.com

Hi Rutger,

[...]

We use these geotiff against a mapserver (and other visualisation
software), and it's working so far. But we use a polar stereographic
projection...

I tried running the same script with only a different area_def, i used the "euro_north" which is provided by default in de areas.def file from mpop. This is a Stereographic projection compared to Lambert equal area for the EASE-grid. Using this Stereographic projection the output Geotiff contains proper metadata and is correctly recognized by other software. So it might my mistake by slightly editing the original EASE-grid parameters or maybe the Lambert projection etc. Overall there is no problem i guess.

It might not be so hard to add an option to the geotiff save function
to support floating values... let's see. Try with the attached file,
and img.save("myimage.tif", floating_point=True) for example.

That works very well, thanks!

The alpha band can be disabled by providing a fill_value that is not
None. Try putting fill_value to (0.0, ) for a grayscale image.

This solves it indeed, i overlooked the exact working of the fill_value at the mpop site.

No, not abusing... Floating point is supported by geotiff, so there's
no reason it shouldn't be supported by geo_image :)

No problem, thanks for you feedback, it's highly valuable to us !
Best regards,
Martin

Since it all works pretty much how i initially intended, i have attached the outcome of both the visual matplotlib output and the geotiff with the floats. The result is very satisfying, thanks.

Rutger

Op dinsdag 6 maart 2012 11:20:42 UTC+1 schreef Martin Raspaud het volgende:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2009, 2011, 2012.
# SMHI,
# Folkborgsvägen 1,
# Norrköping,
# Sweden
# Author(s):

# This file is part of mpop.
# mpop is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# mpop is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with mpop. If not, see <http://www.gnu.org/licenses/>.
"""Module for geographic images.
"""
import os
import Image as pil
import numpy as np
import mpop.imageo.image
from mpop import CONFIG_PATH
from mpop.imageo.logger import LOG
from mpop.utils import ensure_dir
class GeoImage(mpop.imageo.image.Image):
"""This class defines geographic images. As such, it contains not only data
of the different *channels* of the image, but also the area on which it is
defined (*area* parameter) and *time_slot* of the snapshot.

The channels are considered to contain floating point values in the range
[0.0,1.0]. In order to normalize the input data, the *crange* parameter
defines the original range of the data. The conversion to the classical
[0,255] range and byte type is done automagically when saving the image to
file.
See also :class:`image.Image` for more information.
"""
def __init__(self, channels, area, time_slot,
mode = "L", crange = None, fill_value = None, palette = None):
self.area = area
self.time_slot = time_slot
self.tags = {}
self.gdal_options = {}
super(GeoImage, self).__init__(channels, mode, crange,
fill_value, palette)
def save(self, filename, compression=6,
tags=None, gdal_options=None,
fformat=None, blocksize=256, **kwargs):
"""Save the image to the given *filename*. If the extension is "tif",
the image is saved to geotiff_ format, in which case the *compression*
level can be given ([0, 9], 0 meaning off). See also
:meth:`image.Image.save`, :meth:`image.Image.double_save`, and
:meth:`image.Image.secure_save`. The *tags* argument is a dict of tags
to include in the image (as metadata), and the *gdal_options* holds
options for the gdal saving driver. A *blocksize* other than 0 will
result in a tiled image (if possible), with tiles of size equal to
*blocksize*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
file_tuple = os.path.splitext(filename)
fformat = fformat or file_tuple[1][1:]
if fformat.lower() in ('tif', 'tiff'):
self.geotiff_save(filename, compression, tags,
gdal_options, blocksize, **kwargs)
else:
super(GeoImage, self).save(filename, compression, format=fformat, **kwargs)
def _gdal_write_channels(self, dst_ds, channels, opacity, fill_value):
"""Write *channels* in a gdal raster structure *dts_ds*, using
*opacity* as alpha value for valid data, and *fill_value*.
"""
if fill_value is not None:
for i in range(len(channels)):
chn = channels[i].filled(fill_value[i])
dst_ds.GetRasterBand(i + 1).WriteArray(chn)
else:
mask = np.zeros(channels[0].shape, dtype=np.uint8)
i = 0
for i in range(len(channels)):
dst_ds.GetRasterBand(i + 1).WriteArray(channels[i].filled(i))
mask |= np.ma.getmaskarray(channels[i])

try:
mask |= np.ma.getmaskarray(opacity)
except AttributeError:
pass

alpha = np.where(mask, 0, opacity).astype(np.uint8)
dst_ds.GetRasterBand(i + 2).WriteArray(alpha)
def geotiff_save(self, filename, compression=6,
tags=None, gdal_options=None,
blocksize=0, geotransform=None,
spatialref=None, floating_point=False):
"""Save the image to the given *filename* in geotiff_ format, with the
*compression* level in [0, 9]. 0 means not compressed. The *tags*
argument is a dict of tags to include in the image (as metadata). By
default it uses the 'area' instance to generate geotransform and
spatialref information, this can be overwritten by the arguments
*geotransform* and *spatialref*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
from osgeo import gdal, osr

raster = gdal.GetDriverByName("GTiff")
if floating_point:
if self.mode != "L":
raise ValueError("Image must be in 'L' mode for floating point"
" geotif saving")
channels = [self.channels[0].astype(np.float64)]
fill_value = self.fill_value or 0
gformat = gdal.GDT_Float64
else:
channels, fill_value = self._finalize()
gformat = gdal.GDT_Byte
LOG.debug("Saving to GeoTiff.")
if tags is not None:
self.tags.update(tags)
if gdal_options is not None:
self.gdal_options.update(gdal_options)
g_opts = ["=".join(i) for i in self.gdal_options.items()]
if compression != 0:
g_opts.append("COMPRESS=DEFLATE")
g_opts.append("ZLEVEL=" + str(compression))
if blocksize != 0:
g_opts.append("TILED=YES")
g_opts.append("BLOCKXSIZE=" + str(blocksize))
g_opts.append("BLOCKYSIZE=" + str(blocksize))

if(self.mode == "L"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
1,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "LA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds,
channels[:-1], channels[1],
fill_value)
elif(self.mode == "RGB"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
3,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "RGBA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels[:-1], channels[3], fill_value)
else:
raise NotImplementedError("Saving to GeoTIFF using image mode"
" %s is not implemented."%self.mode)

# Create raster GeoTransform based on upper left corner and pixel
# resolution ... if not overwritten by argument geotranform.
if geotransform:
dst_ds.SetGeoTransform(geotransform)
if spatialref:
if not isinstance(spatialref, str):
spatialref = spatialref.ExportToWkt()
dst_ds.SetProjection(spatialref)
else:
try:
from pyresample import utils
from mpop.projector import get_area_def

area = get_area_def(self.area)
except (utils.AreaNotFound, AttributeError):
area = self.area

try:
adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
area.area_extent[3], 0, -area.pixel_size_y]
dst_ds.SetGeoTransform(adfgeotransform)
srs = osr.SpatialReference()
srs.SetProjCS(area.proj_id)
srs.ImportFromProj4(area.proj4_string)
srs = srs.ExportToWkt()
dst_ds.SetProjection(srs)
except AttributeError:
LOG.exception("Could not load geographic data, invalid area")
self.tags.update({'TIFFTAG_DATETIME':
self.time_slot.strftime("%Y:%m:%d %H:%M:%S")})
dst_ds.SetMetadata(self.tags, '')

# Close the dataset

dst_ds = None

def add_overlay(self, color=(0, 0, 0), width=0.5, resolution=None):
"""Add coastline and political borders to image, using *color*.
Loses the masks !

*resolution* is chosen automatically if None (default), otherwise it should be one of:
+-----+-------------------------+---------+
| 'f' | Full resolution | 0.04 km |
|'h' | High resolution | 0.2 km |
|'i' | Intermediate resolution | 1.0 km |
|'l' | Low resolution | 5.0 km |
|'c' | Crude resolution | 25 km |
+-----+-------------------------+---------+
"""

img = self.pil_image()
import ConfigParser
conf = ConfigParser.ConfigParser()
conf.read(os.path.join(CONFIG_PATH, "mpop.cfg"))
coast_dir = conf.get('shapes', 'dir')
LOG.debug("Getting area for overlay: " + str(self.area))
if self.area is None:
raise ValueError("Area of image is None, can't add overlay.")
from mpop.projector import get_area_def
if isinstance(self.area, str):
self.area = get_area_def(self.area)
LOG.info("Add coastlines and political borders to image.")
LOG.debug("Area = " + str(self.area))
if resolution is None:

x_resolution = ((self.area.area_extent[2] -
self.area.area_extent[0]) /
self.area.x_size)
y_resolution = ((self.area.area_extent[3] -
self.area.area_extent[1]) /
self.area.y_size)
res = min(x_resolution, y_resolution)
if res > 25000:
resolution = "c"
elif res > 5000:
resolution = "l"
elif res > 1000:
resolution = "i"
elif res > 200:
resolution = "h"
else:
resolution = "f"
LOG.debug("Automagically choose resolution " + resolution)

from pycoast import ContourWriterAGG
cw_ = ContourWriterAGG(coast_dir)
cw_.add_coastlines(img, self.area, outline=color,
resolution=resolution, width=width)
cw_.add_borders(img, self.area, outline=color,
resolution=resolution, width=width)
arr = np.array(img)
for idx in range(len(self.channels)):
self.channels[idx] = np.ma.array(arr[:, :, idx] / 255.0)

Op dinsdag 6 maart 2012 11:20:42 UTC+1 schreef Martin Raspaud het volgende:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2009, 2011, 2012.
# SMHI,
# Folkborgsvägen 1,
# Norrköping,
# Sweden
# Author(s):
# This file is part of mpop.
# mpop is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# mpop is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with mpop. If not, see <http://www.gnu.org/licenses/>.
"""Module for geographic images.
"""
import os
import Image as pil
import numpy as np
import mpop.imageo.image
from mpop import CONFIG_PATH
from mpop.imageo.logger import LOG
from mpop.utils import ensure_dir
class GeoImage(mpop.imageo.image.Image):
"""This class defines geographic images. As such, it contains not only data
of the different *channels* of the image, but also the area on which it is
defined (*area* parameter) and *time_slot* of the snapshot.

The channels are considered to contain floating point values in the range
[0.0,1.0]. In order to normalize the input data, the *crange* parameter
defines the original range of the data. The conversion to the classical
[0,255] range and byte type is done automagically when saving the image to
file.
See also :class:`image.Image` for more information.
"""
def __init__(self, channels, area, time_slot,
mode = "L", crange = None, fill_value = None, palette = None):
self.area = area
self.time_slot = time_slot
self.tags = {}
self.gdal_options = {}
super(GeoImage, self).__init__(channels, mode, crange,
fill_value, palette)
def save(self, filename, compression=6,
tags=None, gdal_options=None,
fformat=None, blocksize=256, **kwargs):
"""Save the image to the given *filename*. If the extension is "tif",
the image is saved to geotiff_ format, in which case the *compression*
level can be given ([0, 9], 0 meaning off). See also
:meth:`image.Image.save`, :meth:`image.Image.double_save`, and
:meth:`image.Image.secure_save`. The *tags* argument is a dict of tags
to include in the image (as metadata), and the *gdal_options* holds
options for the gdal saving driver. A *blocksize* other than 0 will
result in a tiled image (if possible), with tiles of size equal to
*blocksize*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
file_tuple = os.path.splitext(filename)
fformat = fformat or file_tuple[1][1:]
if fformat.lower() in ('tif', 'tiff'):
self.geotiff_save(filename, compression, tags,
gdal_options, blocksize, **kwargs)
else:
super(GeoImage, self).save(filename, compression, format=fformat, **kwargs)
def _gdal_write_channels(self, dst_ds, channels, opacity, fill_value):
"""Write *channels* in a gdal raster structure *dts_ds*, using
*opacity* as alpha value for valid data, and *fill_value*.
"""
if fill_value is not None:
for i in range(len(channels)):
chn = channels[i].filled(fill_value[i])
dst_ds.GetRasterBand(i + 1).WriteArray(chn)
else:
mask = np.zeros(channels[0].shape, dtype=np.uint8)
i = 0
for i in range(len(channels)):
dst_ds.GetRasterBand(i + 1).WriteArray(channels[i].filled(i))
mask |= np.ma.getmaskarray(channels[i])

try:
mask |= np.ma.getmaskarray(opacity)
except AttributeError:
pass

alpha = np.where(mask, 0, opacity).astype(np.uint8)
dst_ds.GetRasterBand(i + 2).WriteArray(alpha)
def geotiff_save(self, filename, compression=6,
tags=None, gdal_options=None,
blocksize=0, geotransform=None,
spatialref=None, floating_point=False):
"""Save the image to the given *filename* in geotiff_ format, with the
*compression* level in [0, 9]. 0 means not compressed. The *tags*
argument is a dict of tags to include in the image (as metadata). By
default it uses the 'area' instance to generate geotransform and
spatialref information, this can be overwritten by the arguments
*geotransform* and *spatialref*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
from osgeo import gdal, osr

raster = gdal.GetDriverByName("GTiff")
if floating_point:
if self.mode != "L":
raise ValueError("Image must be in 'L' mode for floating point"
" geotif saving")
channels = [self.channels[0].astype(np.float64)]
fill_value = self.fill_value or 0
gformat = gdal.GDT_Float64
else:
channels, fill_value = self._finalize()
gformat = gdal.GDT_Byte
LOG.debug("Saving to GeoTiff.")
if tags is not None:
self.tags.update(tags)
if gdal_options is not None:
self.gdal_options.update(gdal_options)
g_opts = ["=".join(i) for i in self.gdal_options.items()]
if compression != 0:
g_opts.append("COMPRESS=DEFLATE")
g_opts.append("ZLEVEL=" + str(compression))
if blocksize != 0:
g_opts.append("TILED=YES")
g_opts.append("BLOCKXSIZE=" + str(blocksize))
g_opts.append("BLOCKYSIZE=" + str(blocksize))

if(self.mode == "L"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
1,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "LA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds,
channels[:-1], channels[1],
fill_value)
elif(self.mode == "RGB"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
3,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "RGBA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels[:-1], channels[3], fill_value)
else:
raise NotImplementedError("Saving to GeoTIFF using image mode"
" %s is not implemented."%self.mode)

# Create raster GeoTransform based on upper left corner and pixel
# resolution ... if not overwritten by argument geotranform.
if geotransform:
dst_ds.SetGeoTransform(geotransform)
if spatialref:
if not isinstance(spatialref, str):
spatialref = spatialref.ExportToWkt()
dst_ds.SetProjection(spatialref)
else:
try:
from pyresample import utils
from mpop.projector import get_area_def

area = get_area_def(self.area)
except (utils.AreaNotFound, AttributeError):
area = self.area

try:
adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
area.area_extent[3], 0, -area.pixel_size_y]
dst_ds.SetGeoTransform(adfgeotransform)
srs = osr.SpatialReference()
srs.SetProjCS(area.proj_id)
srs.ImportFromProj4(area.proj4_string)
srs = srs.ExportToWkt()
dst_ds.SetProjection(srs)
except AttributeError:
LOG.exception("Could not load geographic data, invalid area")
self.tags.update({'TIFFTAG_DATETIME':
self.time_slot.strftime("%Y:%m:%d %H:%M:%S")})
dst_ds.SetMetadata(self.tags, '')

# Close the dataset

dst_ds = None

def add_overlay(self, color=(0, 0, 0), width=0.5, resolution=None):
"""Add coastline and political borders to image, using *color*.
Loses the masks !

*resolution* is chosen automatically if None (default), otherwise it should be one of:
+-----+-------------------------+---------+
| 'f' | Full resolution | 0.04 km |
|'h' | High resolution | 0.2 km |
|'i' | Intermediate resolution | 1.0 km |
|'l' | Low resolution | 5.0 km |
|'c' | Crude resolution | 25 km |
+-----+-------------------------+---------+
"""

img = self.pil_image()
import ConfigParser
conf = ConfigParser.ConfigParser()
conf.read(os.path.join(CONFIG_PATH, "mpop.cfg"))
coast_dir = conf.get('shapes', 'dir')
LOG.debug("Getting area for overlay: " + str(self.area))
if self.area is None:
raise ValueError("Area of image is None, can't add overlay.")
from mpop.projector import get_area_def
if isinstance(self.area, str):
self.area = get_area_def(self.area)
LOG.info("Add coastlines and political borders to image.")
LOG.debug("Area = " + str(self.area))
if resolution is None:

x_resolution = ((self.area.area_extent[2] -
self.area.area_extent[0]) /
self.area.x_size)
y_resolution = ((self.area.area_extent[3] -
self.area.area_extent[1]) /
self.area.y_size)
res = min(x_resolution, y_resolution)
if res > 25000:
resolution = "c"
elif res > 5000:
resolution = "l"
elif res > 1000:
resolution = "i"
elif res > 200:
resolution = "h"
else:
resolution = "f"
LOG.debug("Automagically choose resolution " + resolution)

from pycoast import ContourWriterAGG
cw_ = ContourWriterAGG(coast_dir)
cw_.add_coastlines(img, self.area, outline=color,
resolution=resolution, width=width)
cw_.add_borders(img, self.area, outline=color,
resolution=resolution, width=width)
arr = np.array(img)
for idx in range(len(self.channels)):
self.channels[idx] = np.ma.array(arr[:, :, idx] / 255.0)

Op dinsdag 6 maart 2012 11:20:42 UTC+1 schreef Martin Raspaud het volgende:

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2009, 2011, 2012.
# SMHI,
# Folkborgsvägen 1,
# Norrköping,
# Sweden
# Author(s):

# Martin Raspaud <martin....@smhi.se>
# Adam Dybbroe <adam.d...@smhi.se>
# Esben S. Nielsen <e...@dmi.dk>
# This file is part of mpop.
# mpop is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# mpop is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
# You should have received a copy of the GNU General Public License
# along with mpop. If not, see <http://www.gnu.org/licenses/>.
"""Module for geographic images.
"""
import os
import Image as pil
import numpy as np
import mpop.imageo.image
from mpop import CONFIG_PATH
from mpop.imageo.logger import LOG
from mpop.utils import ensure_dir
class GeoImage(mpop.imageo.image.Image):
"""This class defines geographic images. As such, it contains not only data
of the different *channels* of the image, but also the area on which it is
defined (*area* parameter) and *time_slot* of the snapshot.

The channels are considered to contain floating point values in the range
[0.0,1.0]. In order to normalize the input data, the *crange* parameter
defines the original range of the data. The conversion to the classical
[0,255] range and byte type is done automagically when saving the image to
file.
See also :class:`image.Image` for more information.
"""
def __init__(self, channels, area, time_slot,
mode = "L", crange = None, fill_value = None, palette = None):
self.area = area
self.time_slot = time_slot
self.tags = {}
self.gdal_options = {}
super(GeoImage, self).__init__(channels, mode, crange,
fill_value, palette)
def save(self, filename, compression=6,
tags=None, gdal_options=None,
fformat=None, blocksize=256, **kwargs):
"""Save the image to the given *filename*. If the extension is "tif",
the image is saved to geotiff_ format, in which case the *compression*
level can be given ([0, 9], 0 meaning off). See also
:meth:`image.Image.save`, :meth:`image.Image.double_save`, and
:meth:`image.Image.secure_save`. The *tags* argument is a dict of tags
to include in the image (as metadata), and the *gdal_options* holds
options for the gdal saving driver. A *blocksize* other than 0 will
result in a tiled image (if possible), with tiles of size equal to
*blocksize*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
file_tuple = os.path.splitext(filename)
fformat = fformat or file_tuple[1][1:]
if fformat.lower() in ('tif', 'tiff'):
self.geotiff_save(filename, compression, tags,
gdal_options, blocksize, **kwargs)
else:
super(GeoImage, self).save(filename, compression, format=fformat, **kwargs)
def _gdal_write_channels(self, dst_ds, channels, opacity, fill_value):
"""Write *channels* in a gdal raster structure *dts_ds*, using
*opacity* as alpha value for valid data, and *fill_value*.
"""
if fill_value is not None:
for i in range(len(channels)):
chn = channels[i].filled(fill_value[i])
dst_ds.GetRasterBand(i + 1).WriteArray(chn)
else:
mask = np.zeros(channels[0].shape, dtype=np.uint8)
i = 0
for i in range(len(channels)):
dst_ds.GetRasterBand(i + 1).WriteArray(channels[i].filled(i))
mask |= np.ma.getmaskarray(channels[i])

try:
mask |= np.ma.getmaskarray(opacity)
except AttributeError:
pass

alpha = np.where(mask, 0, opacity).astype(np.uint8)
dst_ds.GetRasterBand(i + 2).WriteArray(alpha)
def geotiff_save(self, filename, compression=6,
tags=None, gdal_options=None,
blocksize=0, geotransform=None,
spatialref=None, floating_point=False):
"""Save the image to the given *filename* in geotiff_ format, with the
*compression* level in [0, 9]. 0 means not compressed. The *tags*
argument is a dict of tags to include in the image (as metadata). By
default it uses the 'area' instance to generate geotransform and
spatialref information, this can be overwritten by the arguments
*geotransform* and *spatialref*.

.. _geotiff: http://trac.osgeo.org/geotiff/
"""
from osgeo import gdal, osr

raster = gdal.GetDriverByName("GTiff")
if floating_point:
if self.mode != "L":
raise ValueError("Image must be in 'L' mode for floating point"
" geotif saving")
channels = [self.channels[0].astype(np.float64)]
fill_value = self.fill_value or 0
gformat = gdal.GDT_Float64
else:
channels, fill_value = self._finalize()
gformat = gdal.GDT_Byte
LOG.debug("Saving to GeoTiff.")
if tags is not None:
self.tags.update(tags)
if gdal_options is not None:
self.gdal_options.update(gdal_options)
g_opts = ["=".join(i) for i in self.gdal_options.items()]
if compression != 0:
g_opts.append("COMPRESS=DEFLATE")
g_opts.append("ZLEVEL=" + str(compression))
if blocksize != 0:
g_opts.append("TILED=YES")
g_opts.append("BLOCKXSIZE=" + str(blocksize))
g_opts.append("BLOCKYSIZE=" + str(blocksize))

if(self.mode == "L"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
1,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "LA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
2,
gformat,
g_opts)
self._gdal_write_channels(dst_ds,
channels[:-1], channels[1],
fill_value)
elif(self.mode == "RGB"):
ensure_dir(filename)
if fill_value is not None:
dst_ds = raster.Create(filename,
self.width,
self.height,
3,
gformat,
g_opts)
else:
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels, 255, fill_value)
elif(self.mode == "RGBA"):
ensure_dir(filename)
g_opts.append("ALPHA=YES")
dst_ds = raster.Create(filename,
self.width,
self.height,
4,
gformat,
g_opts)
self._gdal_write_channels(dst_ds, channels[:-1], channels[3], fill_value)
else:
raise NotImplementedError("Saving to GeoTIFF using image mode"
" %s is not implemented."%self.mode)

# Create raster GeoTransform based on upper left corner and pixel
# resolution ... if not overwritten by argument geotranform.
if geotransform:
dst_ds.SetGeoTransform(geotransform)
if spatialref:
if not isinstance(spatialref, str):
spatialref = spatialref.ExportToWkt()
dst_ds.SetProjection(spatialref)
else:
try:
from pyresample import utils
from mpop.projector import get_area_def

area = get_area_def(self.area)
except (utils.AreaNotFound, AttributeError):
area = self.area

try:
adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
area.area_extent[3], 0, -area.pixel_size_y]
dst_ds.SetGeoTransform(adfgeotransform)
srs = osr.SpatialReference()
srs.SetProjCS(area.proj_id)
srs.ImportFromProj4(area.proj4_string)
srs = srs.ExportToWkt()
dst_ds.SetProjection(srs)
except AttributeError:
LOG.exception("Could not load geographic data, invalid area")
self.tags.update({'TIFFTAG_DATETIME':
self.time_slot.strftime("%Y:%m:%d %H:%M:%S")})
dst_ds.SetMetadata(self.tags, '')

# Close the dataset

dst_ds = None

def add_overlay(self, color=(0, 0, 0), width=0.5, resolution=None):
"""Add coastline and political borders to image, using *color*.
Loses the masks !

*resolution* is chosen automatically if None (default), otherwise it should be one of:
+-----+-------------------------+---------+
| 'f' | Full resolution | 0.04 km |
|'h' | High resolution | 0.2 km |
|'i' | Intermediate resolution | 1.0 km |
|'l' | Low resolution | 5.0 km |
|'c' | Crude resolution | 25 km |
+-----+-------------------------+---------+
"""

img = self.pil_image()
import ConfigParser
conf = ConfigParser.ConfigParser()
conf.read(os.path.join(CONFIG_PATH, "mpop.cfg"))
coast_dir = conf.get('shapes', 'dir')
LOG.debug("Getting area for overlay: " + str(self.area))
if self.area is None:
raise ValueError("Area of image is None, can't add overlay.")
from mpop.projector import get_area_def
if isinstance(self.area, str):
self.area = get_area_def(self.area)
LOG.info("Add coastlines and political borders to image.")
LOG.debug("Area = " + str(self.area))
if resolution is None:

x_resolution = ((self.area.area_extent[2] -
self.area.area_extent[0]) /
self.area.x_size)
y_resolution = ((self.area.area_extent[3] -
self.area.area_extent[1]) /
self.area.y_size)
res = min(x_resolution, y_resolution)
if res > 25000:
resolution = "c"
elif res > 5000:
resolution = "l"
elif res > 1000:
resolution = "i"
elif res > 200:
resolution = "h"
else:
resolution = "f"
LOG.debug("Automagically choose resolution " + resolution)

from pycoast import ContourWriterAGG
cw_ = ContourWriterAGG(coast_dir)
cw_.add_coastlines(img, self.area, outline=color,
resolution=resolution, width=width)
cw_.add_borders(img, self.area, outline=color,
resolution=resolution, width=width)
arr = np.array(img)
for idx in range(len(self.channels)):
self.channels[idx] = np.ma.array(arr[:, :, idx] / 255.0)

pyresample_output.jpg

Reply all

Reply to author

Forward