How to save image files in geotiff format?

[張伯誠]

Hi, all.

I don't know why there are always some error message in saving a image file in geotiff format.

Here are my error messages and codes:


Warning 6: 'DEFLATE' is an unexpected value for COMPRESS creation option of type string-select.
Warning 6: Driver GTiff does not support ZLEVEL creation option
ERROR 1: Cannot create TIFF file due to missing codec for DEFLATE.

Traceback (most recent call last):
File "npp_image.py", line 97, in <module>
    img.save('/root/npp/pictures/dnb_proj_B3.tif')
 
File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 91, in save

    gdal_options, blocksize, **kwargs)

File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 183, in geotiff_save

    self._gdal_write_channels(dst_ds, channels, 255, fill_value)

File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 102, in _gdal_write_channels

    dst_ds.GetRasterBand(i + 1).WriteArray(chn)

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'


***my code:****


areaid='jap_test'
area_def=get_area_def(areaid)



# plot picture


tslots=[datetime(2013,7,1,4,1),datetime(2013,7,1,4,2),datetime(2013,7,1,4,3),datetime(2013,7,1,4,5),datetime(2013,7,1,4,6),\
datetime(2013,7,1,4,8),datetime(2013,7,1,4,9),datetime(2013,7,1,4,11),datetime(2013,7,1,4,12)]


orbit="00001

global_data=[]

for time_slot in tslots:

    global_data.append(PolarFactory.create_scene("npp","","viirs",time_slot,orbit))


#loading the channels

gldata=[]
for glbd in global_data:

        glbd.load(glbd.image.dnb_rgb.prerequisites|glbd.image.dnb.prerequisites)

        gldata.append(glbd)



glb_data=mpop.scene.assemble_segments(gldata)

# map projection

local_data=glb_data.project(areaid,mode="nearest")


img=local_data.image.dnb()
img.show()
img.save('/root/npp/pictures/dnb_proj.tif')




How could I fix this problem?

Thanks in advance.

BoCheng

PS. I have installed the osgeo module, it's ok to import osgeo.

[Lars Orum Rasmussen]

Hi BoCheng,

Just for isolating the problem, do you have any problems in saving it without any compression ?, like:

img.save('/root/npp/pictures/dnb_proj.tif', compression=0)

    Regards, Lars

2013/7/22 張伯誠 <bomak...@gmail.com>

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[張伯誠]

Hi, Lars.

Thanks for your reply.

I try your code:
img.save('/root/npp/pictures/dnb_proj.tif', compression=0)

There are still some error message as below:

[ERROR: 2013-07-23 14:51:07 : imageo] Could not load geographic data, invalid area

Traceback (most recent call last):

  File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 253, in geotiff_save

    adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,

AttributeError: 'NoneType' object has no attribute 'area_extent'


However, the program did generate the image files in geotiff format!

I'm glad to see this result, after all, there are the image files I want.

But why? Where is this problem? Lacking some modules? Wrong settings/arguments? Others?

Thanks again.

BoCheng

 

 



Lars於 2013年7月22日星期一UTC+8下午6時19分44秒寫道：

[Adam Dybbroe]

BoCheng

I just tried you code example here and it works fine for me, no error messages.

Could you just try save the image as a png?
Do you still get error messages?

I suppose it is an issue with gdal. mpop is using gdal to store images in tiff.
What version of gdal do you have? You can check it like this:

>>> import osgeo
>>> print osgeo.__version__

I have 1.7.3 on RedHat ES 6 and these two rpms:

gdal-python-1.7.3-15.el6.x86_64
gdal-1.7.3-15.el6.x86_64


The code I was using here (using my own data and own area of course) was this:

from mpop.satellites import PolarFactory
import mpop.scene
from datetime import datetime, timedelta
from mpop.projector import get_area_def

tdelta = timedelta(seconds=60)

areaid = 'scan1'
area_def = get_area_def(areaid)

starttime = datetime(2012, 12, 26, 1, 6)
endtime = datetime(2012, 12, 26, 1, 20)
tslot = starttime
orbit="06025"
global_data = []
while tslot < endtime:
    global_data.append(PolarFactory.create_scene("npp", "",
                                                 "viirs", tslot, orbit))
    tslot = tslot + tdelta

# loading the channels

gldata = []
for glbd in global_data:
    try:
        glbd.load(glbd.image.dnb_rgb.prerequisites|glbd.image.dnb.prerequisites)
    except IOError:
        print("No file for that time: %r" % glbd.time_slot)

    gldata.append(glbd)

glb_data = mpop.scene.assemble_segments(gldata)

# map projection
local_data = glb_data.project(areaid, mode="nearest")


img = local_data.image.dnb()
img.show()
img.save('/var/tmp/dnb_proj.tif')



Best regards
Adam

-- 
Adam Dybbroe,
Satellite Remote Sensing Scientist,
Numerical models and Remote Sensing,
Core Services, Swedish Meteorological and Hydrological Institute (SMHI)
www.pytroll.org
nwcsaf.smhi.se
www.smhi.se

[張伯誠]

Adam於 2013年7月23日星期二UTC+8下午7時07分25秒寫道：

BoCheng

I just tried you code example here and it works fine for me, no error messages.

Could you just try save the image as a png?
Do you still get error messages?

Hello, Adam.

It's no any problem to save the image as a png file.
 

I suppose it is an issue with gdal. mpop is using gdal to store images in tiff.
What version of gdal do you have? You can check it like this:

>>> import osgeo
>>> print osgeo.__version__

The version is: 1.9.1
 

I have 1.7.3 on RedHat ES 6 and these two rpms:

gdal-python-1.7.3-15.el6.x86_64
gdal-1.7.3-15.el6.x86_64

I have the module:GDAL-1.9.1-py2.7-linux_x86_64.egg

I guess it may be the version is too new to use?

[Adam Dybbroe]

Hi,

On 2013-07-24 03:59, 張伯誠 wrote:

I guess it may be the version is too new to use?

Yes, maybe there is an issue with mpop and the new version of gdal. We have not yet tested it against your version.
Sorry I didn't pay so much attention to the exact error message, it is quite helpful, see below:

On 2013-07-23 09:25, 張伯誠 wrote:

[ERROR: 2013-07-23 14:51:07 : imageo] Could not load geographic data, invalid area

Traceback (most recent call last):

  File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 253, in geotiff_save

    adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,

AttributeError: 'NoneType' object has no attribute 'area_extent'

Apparently your area is "None" inside geo_image.py

Before you save your image could you just check your area attribute:
print img.area
img.save(...

In addition you might want to add some print statements inside geo_image.py just above the line 253. For example:

...
        else:
            try:
                from pyresample import utils
                from mpop.projector import get_area_def
           
                area = get_area_def(self.area)
                print("Try")
            except (utils.AreaNotFound, AttributeError):
                print("Except")
                area = self.area

            print("self.area: " + str(dir(self.area)))
            print("area: " + str(dir(area)))

            try:

                adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,

                                   area.area_extent[3], 0, -area.pixel_size_y]

...




Also instead of

img.save('/var/tmp/dnb_proj.tif')

Try this:

geotrf = [area_def.area_extent[0], area_def.pixel_size_x, 0,
             area_def.area_extent[3], 0, -area_def.pixel_size_y]
img.save('/var/tmp/dnb_proj.tif', geotransform=geotrf)

Does that work better?


Regards
Adam

[Adam Dybbroe]

BoCheng,

If you follow the alternative path of preparing the tif-tags from
outside mpop, you probably also want to set the projection related tags,
and then do it like this, see mpop and gdal documentation:

geotrf = [area_def.area_extent[0], area_def.pixel_size_x, 0,
area_def.area_extent[3], 0, -area_def.pixel_size_y]

from osgeo import gdal, osr
srs = osr.SpatialReference()
srs.ImportFromProj4(area_def.proj4_string)
srs.SetProjCS(area_def.proj_id)
try:
srs.SetWellKnownGeogCS(area_def.proj_dict['ellps'])
except KeyError:
pass
try:
srs.SetAuthority('PROJCS', 'EPSG',
int(area_def.proj_dict['init'].split('epsg:')[1]))
except (KeyError, IndexError):
pass
srs = srs.ExportToWkt()

img.save('/var/tmp/dnb_proj.tif', geotransform=geotrf, spatialref=srs)


BUT, of course mpop should be able to fix this for you. So we still want
to find out why it doesn't do that in your case.

Regards
Adam

[張伯誠]

Adam於 2013年7月24日星期三UTC+8下午3時35分12秒寫道：

Hi,

Hello, Adam.

Sorry for replying so late because I have another business to do, which is writing a program  can auto-generate satellite image files, it spent lots of my time.

By the way, I work in Centeral Weather Bureau in Taiwan.

On 2013-07-24 03:59, 張伯誠 wrote:

I guess it may be the version is too new to use?

Yes, maybe there is an issue with mpop and the new version of gdal. We have not yet tested it against your version.
Sorry I didn't pay so much attention to the exact error message, it is quite helpful, see below:

On 2013-07-23 09:25, 張伯誠 wrote:

[ERROR: 2013-07-23 14:51:07 : imageo] Could not load geographic data, invalid area

Traceback (most recent call last):

  File "/usr/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py", line 253, in geotiff_save

    adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,

AttributeError: 'NoneType' object has no attribute 'area_extent'

Apparently your area is "None" inside geo_image.py

Before you save your image could you just check your area attribute:
print img.area

It's ok to print img.area.
 

That doesn't work...

My code is the same as before except the line img.save('/root/test.tif',compression=0) which is replaced by  your code as img.save('/root/test.

tif', geotransform=geotrf)

Here are my error messages:


Area ID: jap_test
Name: test
Projection ID: test
Projection: {'a': '6378169.0', 'lat_0': '23.26', 'lon_0': '140', 'proj':
'laea'}
Number of columns: 6800
Number of rows: 4800
Area extent: (-3910243.906307377, -2235711.9561354984, 1707807.495946855,
4236982.8661260195)

Warning 6: 'DEFLATE' is an unexpected value for COMPRESS creation option of
type string-select.
Warning 6: Driver GTiff does not support ZLEVEL creation option
ERROR 1: Cannot create TIFF file due to missing codec for DEFLATE.

Traceback (most recent call last):

  File "make_image.py", line 109, in <module>
    img.save('/root/test.tif', geotransform=geotrf)
  File
"/home/user1/polar2grid/ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 91, in save
    gdal_options, blocksize, **kwargs)
  File

"/home/user1/polar2grid/ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 183, in geotiff_save
    self._gdal_write_channels(dst_ds, channels, 255, fill_value)
  File

"/home/user1/polar2grid/ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 102, in _gdal_write_channels
    dst_ds.GetRasterBand(i + 1).WriteArray(chn)

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

However, don't be worry about that...

I try supplying a extra argument 'compression=0' to the function of img.save(), that is, img.save('/root/test.tif', geotransform=geotrf,compression=0)
(Thanks to Lars' suggestion!)

It generates the geotiff file without any error messages!

Thank you and Lars again!

BoCheng
 

[張伯誠]

Adam於 2013年7月24日星期三UTC+8下午3時55分40秒寫道：

BoCheng,

Hello, Adam.

If you follow the alternative path of preparing the tif-tags from
outside mpop, you probably also want to set the projection related tags,
and then do it like this, see mpop and gdal documentation:

Is your meaning that adding the below code to geo_image.py  around the line 253?
 

[Adam Dybbroe]

Hi,

No, no!
That code was a proposal to your script and an alternative way of using mpop, specifying the tiff-stuff outside mpop. You just tried adding the geotransform, and you said it failed, so probably this too will fail, but the example above is more complete. Here I get the same result (same tiff file) as if doing img.save('/var/tmp/dnb_proj.tif').

My suggestion in geo_image.py was to add some print statements just to understand what is going on inside geo_image.py in your case. So in order for me to better understand your problem make these changes around line 253:

...
        else:
            try:
                from pyresample import utils
                from mpop.projector import get_area_def
           
                area = get_area_def(self.area)
                print("Try")
            except (utils.AreaNotFound, AttributeError):
                print("Except")
                area = self.area

            print("self.area: " + str(dir(self.area)))
            print("area: " + str(dir(area)))

            try:
                adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
                                   area.area_extent[3], 0, -area.pixel_size_y]

...

and just sent me the output.


Regards
Adam

BUT, of course mpop should be able to fix this for you. So we still want
to find out why it doesn't do that in your case.

Regards
Adam

--
You received this message because you are subscribed to the Google Groups "pytroll" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pytroll+u...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.
 
 

[張伯誠]

Adam於 2013年7月26日星期五UTC+8下午4時24分48秒寫道：

OK, I understand your meaning.
 

My suggestion in geo_image.py was to add some print statements just to understand what is going on inside geo_image.py in your case. So in order for me to better understand your problem make these changes around line 253:

...
        else:
            try:
                from pyresample import utils
                from mpop.projector import get_area_def
           
                area = get_area_def(self.area)
                print("Try")
            except (utils.AreaNotFound, AttributeError):
                print("Except")
                area = self.area

            print("self.area: " + str(dir(self.area)))
            print("area: " + str(dir(area)))

            try:
                adfgeotransform = [area.area_extent[0], area.pixel_size_x, 0,
                                   area.area_extent[3], 0, -area.pixel_size_y]

...

and just sent me the output.

 

 I have added the code above to geo_image.py around line 253.

 but it doesn't print anything, for example, "Try" , "Except", "self.area:", or "area:"


Here are my erroe message by using img.save('/root/test.tif', geotransform=geotrf):

Warning 6: 'DEFLATE' is an unexpected value for COMPRESS creation option of

type string-select.
Warning 6: Driver GTiff does not support ZLEVEL creation option
ERROR 1: Cannot create TIFF file due to missing codec for DEFLATE.

Traceback (most recent call last):

  File "make_image.py", line 109, in <module>
    img.save('/root/test.tif', geotransform=geotrf)
  File
"/home/user1/polar2grid/

ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 91, in save
    gdal_options, blocksize, **kwargs)
  File

"/home/user1/polar2grid/ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 183, in geotiff_save
    self._gdal_write_channels(dst_ds, channels, 255, fill_value)
  File

"/home/user1/polar2grid/ShellB3/lib/python2.7/site-packages/mpop-0.13.1-py2.7.egg/mpop/imageo/geo_image.py",

line 102, in _gdal_write_channels
    dst_ds.GetRasterBand(i + 1).WriteArray(chn)

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'


BoCheng

[張伯誠]

Hi, Adam.

I guess it may be the problem is about the parameter 'compression'.

Because as I add a extra argument 'compression=0' to the img.save(), which is img.save('/root/test.tif',compression=0,geotransform=geotrf) , there are no any error messages in saving files in geotiff format.

Thanks for your help.

BoCheng





張伯誠於 2013年7月26日星期五UTC+8下午4時52分19秒寫道：

[張伯誠]

Hello, everybody:

I have a new problem: can I save a image file with coastlines and lat-lon lines in geotiff format?

Thank Adam that it has no any error messages in saving image files in geotiff format, and the image files do have the informations about parameters of projection, area_def, pixel-size, latitude-longitude, time-tag, and etc.(by the software OpenEV 1.80)

But the image file with coastlines and lat-lon lines and saved in geotiff format has no the informatin metioned above.

Here is my code:

**********************************************

from PIL import Image, ImageFont
from pycoast import ContourWriter
from pycoast import ContourWriterAGG
from mpop.projector import get_area_def
from mpop.imageo import geo_image
from mpop.satellites import PolarFactory
from datetime import datetime,timedelta
import mpop

areaid='jap_test'
area_def=get_area_def(areaid)



tslots=[]
start_time=datetime(2013,7,1,4,1)
tstep=timedelta(seconds=60)
end_time=start_time+timedelta(minutes=20)
orbit="00001"

while start_time<end_time:
    tslots.append(start_time)
    start_time+=tstep

global_data=[]

for time_slot in tslots:
    global_data.append(PolarFactory.create_scene("npp","","viirs",time_slot,orbit))


#loading the channels

gldata=[]
for glbd in global_data:

    try:
        glbd.load(['DNB'])
        #glbd.area=glbd['DNB'].area
            gldata.append(glbd)
    except:
        pass

global_data=mpop.scene.assemble_segments(global_data)

local_data=global_data.project(areaid,mode="nearest")
img=geo_image.GeoImage(local_data['DNB'].data,areaid,tslots[0],fill_value=0,mode='L')
img.enhance(stretch='linear')

geotrf = [area_def.area_extent[0], area_def.pixel_size_x, 0,
           area_def.area_extent[3], 0, -area_def.pixel_size_y]


from osgeo import gdal, osr
srs = osr.SpatialReference()
srs.ImportFromProj4(area_def.proj4_string)
srs.SetProjCS(area_def.proj_id)
try:
     srs.SetWellKnownGeogCS(area_def.proj_dict['ellps'])
except KeyError:
     pass
try:
     srs.SetAuthority('PROJCS', 'EPSG',
int(area_def.proj_dict['init'].split('epsg:')[1]))
except (KeyError, IndexError):
     pass
srs = srs.ExportToWkt()

img.save('/root/test.tif', geotransform=geotrf,spatialref=srs,compression=0)


font=ImageFont.truetype('/usr/share/fonts/dejavu/DejaVuSerif.ttf',64)
cw=ContourWriter('/home/user1/GSHHS_DATA_ROOT')
cw.add_coastlines_to_file('/root/test.tif',area_def,resolution='l',level=4)
cw.add_grid_to_file('/root/test_other.tif',area_def,(10.0,10.0),(2.0,2.0),font,fill='red',outline='red',minor_outline='red')

***********************

I try different ways, the results are the same: missing the informations of geotiff files should have.

1.As above code, save a projected image in geotiff format, and then add coastlines and lat-lon lines.
2.save a projected image in geotiff format, and then open it, adding coastlines and lat-lon lines, saving in geotiff format with different file name.
3.save a projected image in png format, and then open it, adding coastlines and lat-lon lines, saving in geotiff format with different file name.

How should I do it?

Thanks for your help.

Bocheng

[Martin Raspaud]

Hi Bocheng,

> I try different ways, the results are the same: missing the informations of
> geotiff files should have.
>
> 1.As above code, save a projected image in geotiff format, and then add
> coastlines and lat-lon lines.
> 2.save a projected image in geotiff format, and then open it, adding
> coastlines and lat-lon lines, saving in geotiff format with different file
> name.
> 3.save a projected image in png format, and then open it, adding coastlines
> and lat-lon lines, saving in geotiff format with different file name.
>
> How should I do it?

The problem is that pycoast uses PIL to read and write tiff files, and
doesn't know anything about geotif, so I guess it just discards the geo
info from the file when saving it.

I have a small method called that is called add_overlay for the geoimage
objects (geo_image.py file) that works around that. It is probably not
enough for your needs, but look at the code, you'll surely figure out
how to do then.

Best regards,
Martin

[張伯誠]

Hello, Martin.  Long time no see!

I use the method add_overlay(), it plots coastlines on image and keeps the geographic informations! That's what I want!

I'll try adding lat-lon grids by the code at add_overlay().

Thank you, Martin!



Martin Raspaud於 2013年8月1日星期四UTC+8下午3時08分21秒寫道：