[GRASS-user] Import/conversion problem

Adam_Dershowitz · February 11, 2011, 5:23pm

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.
I downloaded some aerial imagery from here:
http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211--,00.html

I want to import it into a UTM zone 17 project. So I did the following:
gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' input.sid output_utm.sid

gdapwarp responds with an error, but then continues with the conversion:
ERROR 6: Failed to initialize PROJ.4 with `+proj=omerc +lat_0=0 +lonc=0 +alpha=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs '.
lat_0 = 0 or 90 or alpha = 90
0...10...20...30...40...50...60...70...80...90...100 - done.

If I then try to import the image, into my project, it does work, but places it WAY off from where it should be (millions of meters!). I guess I should not be surprised, since the error says that it is trying to use a lat/long of 0. But I don't understand why it is not working correctly.

So, clearly something is wrong with how I am doing the conversion. Am I missing a flag or something? It looks to me like gdalinfo is doing a correct read, but somehow gdalwarp is missing something on the input.
Any help would be greatly appreciated.

If I do gdalinfo on the input file, here is what I get:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: oakgrove_sw.sid
       oakgrove_sw.sdw
       oakgrove_sw.sid.aux.xml
Size is 6415, 8061
Coordinate System is:
PROJCS["IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010002,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Hotine_Oblique_Mercator"],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",0],
    PARAMETER["azimuth",0],
    PARAMETER["rectified_grid_angle",90],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (663144.500000000000000,240035.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  IMAGE__INPUT_NAME=K:\blk18\oakgrove_sw.tif
  IMAGE__INPUT_FILE_SIZE=155199578.000000
  GEOTIFF_NUM__1024__GTModelTypeGeoKey=1
  GEOTIFF_CHAR__GTModelTypeGeoKey=ModelTypeProjected
  GEOTIFF_NUM__1025__GTRasterTypeGeoKey=1
  GEOTIFF_CHAR__GTRasterTypeGeoKey=RasterPixelIsArea
  GEOTIFF_NUM__1026__GTCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
  GEOTIFF_NUM__2048__GeographicTypeGeoKey=4269
  GEOTIFF_CHAR__GeographicTypeGeoKey=GCS_NAD83
  GEOTIFF_NUM__2060__GeogAzimuthUnitsGeoKey=9102
  GEOTIFF_CHAR__GeogAzimuthUnitsGeoKey=Angular_Degree
  GEOTIFF_NUM__3072__ProjectedCSTypeGeoKey=32767
  GEOTIFF_CHAR__ProjectedCSTypeGeoKey=User-Defined
  GEOTIFF_NUM__3073__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
  GEOTIFF_NUM__3074__ProjectionGeoKey=32767
  GEOTIFF_CHAR__ProjectionGeoKey=User-Defined
  GEOTIFF_NUM__3075__ProjCoordTransGeoKey=3
  GEOTIFF_CHAR__ProjCoordTransGeoKey=CT_ObliqueMercator
  GEOTIFF_NUM__3076__ProjLinearUnitsGeoKey=9001
  GEOTIFF_CHAR__ProjLinearUnitsGeoKey=Linear_Meter
  GEOTIFF_NUM__3088__ProjCenterLongGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  GEOTIFF_NUM__3089__ProjCenterLatGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  GEOTIFF_NUM__3090__ProjCenterEastingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  GEOTIFF_NUM__3091__ProjCenterNorthingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  GEOTIFF_NUM__3093__ProjScaleAtCenterGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  GEOTIFF_NUM__3094__ProjAzimuthAngleGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
  IMAGE__Z_RESOLUTION=0.000000
  GEO__ModelTypeGeoKey=1
  GEO__ProjectedCSTypeGeoKey=32767
  GEO__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
  IMAGE__Z_ORIGIN=0.000000
  IMAGE__COMPRESSION_VERSION=1,6,1
  IMAGE__TARGET_COMPRESSION_RATIO=20.000000
  IMAGE__COMPRESSION_NLEV=4
  IMAGE__COMPRESSION_WEIGHT=4.000000
  IMAGE__COMPRESSION_GAMMA=2.000000
  IMAGE__COMPRESSION_BLOCK_SIZE=512
  IMAGE__CREATION_DATE=Wed Mar 10 17:11:37 2004

  VERSION=MG2
  IMAGE__BITS_PER_SAMPLE=8
  IMAGE__WIDTH=6415
  IMAGE__HEIGHT=8061
  IMAGE__COLOR_SCHEME=0
  IMAGE__DATA_TYPE=0
  IMAGE__XY_ORIGIN=663145.000000,240035.000000
  IMAGE__X_RESOLUTION=1.000000
  IMAGE__Y_RESOLUTION=1.000000
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 663144.500, 240035.500)
Lower Left ( 663144.500, 231974.500)
Upper Right ( 669559.500, 240035.500)
Lower Right ( 669559.500, 231974.500)
Center ( 666352.000, 236005.000)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Red
  Minimum=0.000, Maximum=255.000, Mean=97.677, StdDev=49.163
  Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 2 Block=1024x128 Type=Byte, ColorInterp=Green
  Minimum=0.000, Maximum=255.000, Mean=97.555, StdDev=49.371
  Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue
  Minimum=0.000, Maximum=255.000, Mean=97.402, StdDev=49.024
  Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

--Adam

Stu1 · February 11, 2011, 6:47pm

MrSID files are a bit tricky. This is a proprietary format owned by LizardTech (see http://www.gdal.org/frmt_mrsid.html). That being said, they freely provide decode tools at their website (http://www.lizardtech.com/). Working on OS X, I use their Raster_DSDK in a little command line routine to convert the file to a geotiff (tifg) and then import it into GRASS with r.in.gdal. Note that the MrSID format is very efficient and boils a quad down to about 1.5 Gb. This will become 3 or 4 times greater as a tif so make sure you have plenty of room to store it and process it. I use an old version of Expressview (a LizardTech viewer that used to be available for OS X) to look at the image and determine the coordinates of the bit that I really want in order to minimize the resulting file sizes.

AFAIK there are no open source apps that deal with MrSID directly (for those that do, you have to install your own version of the decoder first)

Stu

On Feb 11, 2011, at 12:23 PM, Adam Dershowitz, Ph.D., P.E. wrote:

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.
I downloaded some aerial imagery from here:
http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211--,00.html

I want to import it into a UTM zone 17 project. So I did the following:
gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' input.sid output_utm.sid

gdapwarp responds with an error, but then continues with the conversion:
ERROR 6: Failed to initialize PROJ.4 with `+proj=omerc +lat_0=0 +lonc=0 +alpha=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs '.
lat_0 = 0 or 90 or alpha = 90
0...10...20...30...40...50...60...70...80...90...100 - done.

If I then try to import the image, into my project, it does work, but places it WAY off from where it should be (millions of meters!). I guess I should not be surprised, since the error says that it is trying to use a lat/long of 0. But I don't understand why it is not working correctly.

So, clearly something is wrong with how I am doing the conversion. Am I missing a flag or something? It looks to me like gdalinfo is doing a correct read, but somehow gdalwarp is missing something on the input.
Any help would be greatly appreciated.

If I do gdalinfo on the input file, here is what I get:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: oakgrove_sw.sid
      oakgrove_sw.sdw
      oakgrove_sw.sid.aux.xml
Size is 6415, 8061
Coordinate System is:
PROJCS["IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE",
   GEOGCS["NAD83",
       DATUM["North_American_Datum_1983",
           SPHEROID["GRS 1980",6378137,298.2572221010002,
               AUTHORITY["EPSG","7019"]],
           AUTHORITY["EPSG","6269"]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433],
       AUTHORITY["EPSG","4269"]],
   PROJECTION["Hotine_Oblique_Mercator"],
   PARAMETER["latitude_of_center",0],
   PARAMETER["longitude_of_center",0],
   PARAMETER["azimuth",0],
   PARAMETER["rectified_grid_angle",90],
   PARAMETER["scale_factor",1],
   PARAMETER["false_easting",0],
   PARAMETER["false_northing",0],
   UNIT["metre",1,
       AUTHORITY["EPSG","9001"]]]
Origin = (663144.500000000000000,240035.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
IMAGE__INPUT_NAME=K:\blk18\oakgrove_sw.tif
IMAGE__INPUT_FILE_SIZE=155199578.000000
GEOTIFF_NUM__1024__GTModelTypeGeoKey=1
GEOTIFF_CHAR__GTModelTypeGeoKey=ModelTypeProjected
GEOTIFF_NUM__1025__GTRasterTypeGeoKey=1
GEOTIFF_CHAR__GTRasterTypeGeoKey=RasterPixelIsArea
GEOTIFF_NUM__1026__GTCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
GEOTIFF_NUM__2048__GeographicTypeGeoKey=4269
GEOTIFF_CHAR__GeographicTypeGeoKey=GCS_NAD83
GEOTIFF_NUM__2060__GeogAzimuthUnitsGeoKey=9102
GEOTIFF_CHAR__GeogAzimuthUnitsGeoKey=Angular_Degree
GEOTIFF_NUM__3072__ProjectedCSTypeGeoKey=32767
GEOTIFF_CHAR__ProjectedCSTypeGeoKey=User-Defined
GEOTIFF_NUM__3073__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
GEOTIFF_NUM__3074__ProjectionGeoKey=32767
GEOTIFF_CHAR__ProjectionGeoKey=User-Defined
GEOTIFF_NUM__3075__ProjCoordTransGeoKey=3
GEOTIFF_CHAR__ProjCoordTransGeoKey=CT_ObliqueMercator
GEOTIFF_NUM__3076__ProjLinearUnitsGeoKey=9001
GEOTIFF_CHAR__ProjLinearUnitsGeoKey=Linear_Meter
GEOTIFF_NUM__3088__ProjCenterLongGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3089__ProjCenterLatGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3090__ProjCenterEastingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3091__ProjCenterNorthingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3093__ProjScaleAtCenterGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3094__ProjAzimuthAngleGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
IMAGE__Z_RESOLUTION=0.000000
GEO__ModelTypeGeoKey=1
GEO__ProjectedCSTypeGeoKey=32767
GEO__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
IMAGE__Z_ORIGIN=0.000000
IMAGE__COMPRESSION_VERSION=1,6,1
IMAGE__TARGET_COMPRESSION_RATIO=20.000000
IMAGE__COMPRESSION_NLEV=4
IMAGE__COMPRESSION_WEIGHT=4.000000
IMAGE__COMPRESSION_GAMMA=2.000000
IMAGE__COMPRESSION_BLOCK_SIZE=512
IMAGE__CREATION_DATE=Wed Mar 10 17:11:37 2004

VERSION=MG2
IMAGE__BITS_PER_SAMPLE=8
IMAGE__WIDTH=6415
IMAGE__HEIGHT=8061
IMAGE__COLOR_SCHEME=0
IMAGE__DATA_TYPE=0
IMAGE__XY_ORIGIN=663145.000000,240035.000000
IMAGE__X_RESOLUTION=1.000000
IMAGE__Y_RESOLUTION=1.000000
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 663144.500, 240035.500)
Lower Left ( 663144.500, 231974.500)
Upper Right ( 669559.500, 240035.500)
Lower Right ( 669559.500, 231974.500)
Center ( 666352.000, 236005.000)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Red
Minimum=0.000, Maximum=255.000, Mean=97.677, StdDev=49.163
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 2 Block=1024x128 Type=Byte, ColorInterp=Green
Minimum=0.000, Maximum=255.000, Mean=97.555, StdDev=49.371
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue
Minimum=0.000, Maximum=255.000, Mean=97.402, StdDev=49.024
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

--Adam

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Adam_Dershowitz · February 11, 2011, 7:08pm

I am also working on OS X. But, I am using the MrSID plugin, from William Kyngesburye. So it seems that both gdalwarp and r.in.gdal are willing to import it. But, maybe I will try to search around for the apps that you are referring to, to see if they might do a better job on the conversion.
Are you suggesting that the MrSID plugin that I have is working, but is loosing reference data along the way?

--Adam

On Feb 11, 2011, at 10:47 AM, Stuart Edwards wrote:

MrSID files are a bit tricky. This is a proprietary format owned by LizardTech (see http://www.gdal.org/frmt_mrsid.html). That being said, they freely provide decode tools at their website (http://www.lizardtech.com/). Working on OS X, I use their Raster_DSDK in a little command line routine to convert the file to a geotiff (tifg) and then import it into GRASS with r.in.gdal. Note that the MrSID format is very efficient and boils a quad down to about 1.5 Gb. This will become 3 or 4 times greater as a tif so make sure you have plenty of room to store it and process it. I use an old version of Expressview (a LizardTech viewer that used to be available for OS X) to look at the image and determine the coordinates of the bit that I really want in order to minimize the resulting file sizes.

AFAIK there are no open source apps that deal with MrSID directly (for those that do, you have to install your own version of the decoder first)

Stu

On Feb 11, 2011, at 12:23 PM, Adam Dershowitz, Ph.D., P.E. wrote:

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.
I downloaded some aerial imagery from here:
http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211--,00.html

I want to import it into a UTM zone 17 project. So I did the following:
gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' input.sid output_utm.sid

gdapwarp responds with an error, but then continues with the conversion:
ERROR 6: Failed to initialize PROJ.4 with `+proj=omerc +lat_0=0 +lonc=0 +alpha=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs '.
lat_0 = 0 or 90 or alpha = 90
0...10...20...30...40...50...60...70...80...90...100 - done.

If I then try to import the image, into my project, it does work, but places it WAY off from where it should be (millions of meters!). I guess I should not be surprised, since the error says that it is trying to use a lat/long of 0. But I don't understand why it is not working correctly.

So, clearly something is wrong with how I am doing the conversion. Am I missing a flag or something? It looks to me like gdalinfo is doing a correct read, but somehow gdalwarp is missing something on the input.
Any help would be greatly appreciated.

If I do gdalinfo on the input file, here is what I get:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: oakgrove_sw.sid
     oakgrove_sw.sdw
     oakgrove_sw.sid.aux.xml
Size is 6415, 8061
Coordinate System is:
PROJCS["IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE",
  GEOGCS["NAD83",
      DATUM["North_American_Datum_1983",
          SPHEROID["GRS 1980",6378137,298.2572221010002,
              AUTHORITY["EPSG","7019"]],
          AUTHORITY["EPSG","6269"]],
      PRIMEM["Greenwich",0],
      UNIT["degree",0.0174532925199433],
      AUTHORITY["EPSG","4269"]],
  PROJECTION["Hotine_Oblique_Mercator"],
  PARAMETER["latitude_of_center",0],
  PARAMETER["longitude_of_center",0],
  PARAMETER["azimuth",0],
  PARAMETER["rectified_grid_angle",90],
  PARAMETER["scale_factor",1],
  PARAMETER["false_easting",0],
  PARAMETER["false_northing",0],
  UNIT["metre",1,
      AUTHORITY["EPSG","9001"]]]
Origin = (663144.500000000000000,240035.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
IMAGE__INPUT_NAME=K:\blk18\oakgrove_sw.tif
IMAGE__INPUT_FILE_SIZE=155199578.000000
GEOTIFF_NUM__1024__GTModelTypeGeoKey=1
GEOTIFF_CHAR__GTModelTypeGeoKey=ModelTypeProjected
GEOTIFF_NUM__1025__GTRasterTypeGeoKey=1
GEOTIFF_CHAR__GTRasterTypeGeoKey=RasterPixelIsArea
GEOTIFF_NUM__1026__GTCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
GEOTIFF_NUM__2048__GeographicTypeGeoKey=4269
GEOTIFF_CHAR__GeographicTypeGeoKey=GCS_NAD83
GEOTIFF_NUM__2060__GeogAzimuthUnitsGeoKey=9102
GEOTIFF_CHAR__GeogAzimuthUnitsGeoKey=Angular_Degree
GEOTIFF_NUM__3072__ProjectedCSTypeGeoKey=32767
GEOTIFF_CHAR__ProjectedCSTypeGeoKey=User-Defined
GEOTIFF_NUM__3073__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
GEOTIFF_NUM__3074__ProjectionGeoKey=32767
GEOTIFF_CHAR__ProjectionGeoKey=User-Defined
GEOTIFF_NUM__3075__ProjCoordTransGeoKey=3
GEOTIFF_CHAR__ProjCoordTransGeoKey=CT_ObliqueMercator
GEOTIFF_NUM__3076__ProjLinearUnitsGeoKey=9001
GEOTIFF_CHAR__ProjLinearUnitsGeoKey=Linear_Meter
GEOTIFF_NUM__3088__ProjCenterLongGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3089__ProjCenterLatGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3090__ProjCenterEastingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3091__ProjCenterNorthingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3093__ProjScaleAtCenterGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
GEOTIFF_NUM__3094__ProjAzimuthAngleGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000
IMAGE__Z_RESOLUTION=0.000000
GEO__ModelTypeGeoKey=1
GEO__ProjectedCSTypeGeoKey=32767
GEO__PCSCitationGeoKey=IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support
Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $
Projection Name = Oblique Mercator (Hotine)
Units = meters
GeoTIFF Units = meters|
IMAGE__Z_ORIGIN=0.000000
IMAGE__COMPRESSION_VERSION=1,6,1
IMAGE__TARGET_COMPRESSION_RATIO=20.000000
IMAGE__COMPRESSION_NLEV=4
IMAGE__COMPRESSION_WEIGHT=4.000000
IMAGE__COMPRESSION_GAMMA=2.000000
IMAGE__COMPRESSION_BLOCK_SIZE=512
IMAGE__CREATION_DATE=Wed Mar 10 17:11:37 2004

VERSION=MG2
IMAGE__BITS_PER_SAMPLE=8
IMAGE__WIDTH=6415
IMAGE__HEIGHT=8061
IMAGE__COLOR_SCHEME=0
IMAGE__DATA_TYPE=0
IMAGE__XY_ORIGIN=663145.000000,240035.000000
IMAGE__X_RESOLUTION=1.000000
IMAGE__Y_RESOLUTION=1.000000
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 663144.500, 240035.500)
Lower Left ( 663144.500, 231974.500)
Upper Right ( 669559.500, 240035.500)
Lower Right ( 669559.500, 231974.500)
Center ( 666352.000, 236005.000)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Red
Minimum=0.000, Maximum=255.000, Mean=97.677, StdDev=49.163
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 2 Block=1024x128 Type=Byte, ColorInterp=Green
Minimum=0.000, Maximum=255.000, Mean=97.555, StdDev=49.371
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504
Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue
Minimum=0.000, Maximum=255.000, Mean=97.402, StdDev=49.024
Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

--Adam

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Stu1 · February 11, 2011, 8:14pm

I’m not familiar with the Kyngesburye plugin and haven’t used gdalwarp, but here:

I want to import it into a UTM zone 17 project. So I did the following:

gdalwarp -t_srs ‘+proj=utm +zone=17 +datum=WGS84’ input.sid output_utm.sid

it looks like you are just trying to create a new .sid file in zone 17 (output_utm.sid). From the meta file at the MI source, I think the file is already projected in utm - probably in zone 17, and with units of ‘meters’. So what you maybe intended was to just create a .tiffg from the .sid so that it can be imported into GRASS.

In mrsiddecode (the app that does all the work inside the plugin - I assume) I would just issue the command

./mrsiddecode -i Hancock.sid -o Hancock_south.tif -of tifg -ulxy 0 0 -lrxy 71145 29160

where

-i - input
-o - output
-of - output file type
ulxy - upper left x,y
lrxy - lower right x,y

and you get a nice geotiff that r.in.gdal will recognize - already in utm zone 17 most likely. Actually, in my case, the OH mrsids are projected in state plane coordinates but you get the picture.

To get the mrsiddecode app you must register at

http://developer.lizardtech.com/

so you can access to the Mac version of the decoder and mrsidinfo download page. And be slightly comfortable on the command line.

Stu

On Feb 11, 2011, at 2:08 PM, Adam Dershowitz, Ph.D., P.E. wrote:

I am also working on OS X. But, I am using the MrSID plugin, from William Kyngesburye. So it seems that both gdalwarp and r.in.gdal are willing to import it. But, maybe I will try to search around for the apps that you are referring to, to see if they might do a better job on the conversion.
Are you suggesting that the MrSID plugin that I have is working, but is loosing reference data along the way?

–Adam

On Feb 11, 2011, at 10:47 AM, Stuart Edwards wrote:

MrSID files are a bit tricky. This is a proprietary format owned by LizardTech (see http://www.gdal.org/frmt_mrsid.html). That being said, they freely provide decode tools at their website (http://www.lizardtech.com/). Working on OS X, I use their Raster_DSDK in a little command line routine to convert the file to a geotiff (tifg) and then import it into GRASS with r.in.gdal. Note that the MrSID format is very efficient and boils a quad down to about 1.5 Gb. This will become 3 or 4 times greater as a tif so make sure you have plenty of room to store it and process it. I use an old version of Expressview (a LizardTech viewer that used to be available for OS X) to look at the image and determine the coordinates of the bit that I really want in order to minimize the resulting file sizes.

AFAIK there are no open source apps that deal with MrSID directly (for those that do, you have to install your own version of the decoder first)

Stu

On Feb 11, 2011, at 12:23 PM, Adam Dershowitz, Ph.D., P.E. wrote:

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.

I downloaded some aerial imagery from here:

http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211–,00.html

I want to import it into a UTM zone 17 project. So I did the following:

gdalwarp -t_srs ‘+proj=utm +zone=17 +datum=WGS84’ input.sid output_utm.sid

gdapwarp responds with an error, but then continues with the conversion:

ERROR 6: Failed to initialize PROJ.4 with `+proj=omerc +lat_0=0 +lonc=0 +alpha=0 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs '.

lat_0 = 0 or 90 or alpha = 90

0…10…20…30…40…50…60…70…80…90…100 - done.

If I then try to import the image, into my project, it does work, but places it WAY off from where it should be (millions of meters!). I guess I should not be surprised, since the error says that it is trying to use a lat/long of 0. But I don’t understand why it is not working correctly.

So, clearly something is wrong with how I am doing the conversion. Am I missing a flag or something? It looks to me like gdalinfo is doing a correct read, but somehow gdalwarp is missing something on the input.

Any help would be greatly appreciated.

If I do gdalinfo on the input file, here is what I get:

Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)

Files: oakgrove_sw.sid

oakgrove_sw.sdw

oakgrove_sw.sid.aux.xml

Size is 6415, 8061

Coordinate System is:

PROJCS["IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection = Oblique Mercator (Hotine)|IMAGINE",

GEOGCS[“NAD83”,

DATUM[“North_American_Datum_1983”,

SPHEROID[“GRS 1980”,6378137,298.2572221010002,

AUTHORITY[“EPSG”,“7019”]],

AUTHORITY[“EPSG”,“6269”]],

PRIMEM[“Greenwich”,0],

UNIT[“degree”,0.0174532925199433],

AUTHORITY[“EPSG”,“4269”]],

PROJECTION[“Hotine_Oblique_Mercator”],

PARAMETER[“latitude_of_center”,0],

PARAMETER[“longitude_of_center”,0],

PARAMETER[“azimuth”,0],

PARAMETER[“rectified_grid_angle”,90],

PARAMETER[“scale_factor”,1],

PARAMETER[“false_easting”,0],

PARAMETER[“false_northing”,0],

UNIT[“metre”,1,

AUTHORITY[“EPSG”,“9001”]]]

Origin = (663144.500000000000000,240035.500000000000000)

Pixel Size = (1.000000000000000,-1.000000000000000)

Metadata:

IMAGE__INPUT_NAME=K:\blk18\oakgrove_sw.tif

IMAGE__INPUT_FILE_SIZE=155199578.000000

GEOTIFF_NUM__1024__GTModelTypeGeoKey=1

GEOTIFF_CHAR__GTModelTypeGeoKey=ModelTypeProjected

GEOTIFF_NUM__1025__GTRasterTypeGeoKey=1

GEOTIFF_CHAR__GTRasterTypeGeoKey=RasterPixelIsArea

GEOTIFF_NUM__1026__GTCitationGeoKey=IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection Name = Oblique Mercator (Hotine)

Units = meters

GeoTIFF Units = meters|

GEOTIFF_NUM__2048__GeographicTypeGeoKey=4269

GEOTIFF_CHAR__GeographicTypeGeoKey=GCS_NAD83

GEOTIFF_NUM__2060__GeogAzimuthUnitsGeoKey=9102

GEOTIFF_CHAR__GeogAzimuthUnitsGeoKey=Angular_Degree

GEOTIFF_NUM__3072__ProjectedCSTypeGeoKey=32767

GEOTIFF_CHAR__ProjectedCSTypeGeoKey=User-Defined

GEOTIFF_NUM__3073__PCSCitationGeoKey=IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection Name = Oblique Mercator (Hotine)

Units = meters

GeoTIFF Units = meters|

GEOTIFF_NUM__3074__ProjectionGeoKey=32767

GEOTIFF_CHAR__ProjectionGeoKey=User-Defined

GEOTIFF_NUM__3075__ProjCoordTransGeoKey=3

GEOTIFF_CHAR__ProjCoordTransGeoKey=CT_ObliqueMercator

GEOTIFF_NUM__3076__ProjLinearUnitsGeoKey=9001

GEOTIFF_CHAR__ProjLinearUnitsGeoKey=Linear_Meter

GEOTIFF_NUM__3088__ProjCenterLongGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

GEOTIFF_NUM__3089__ProjCenterLatGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

GEOTIFF_NUM__3090__ProjCenterEastingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

GEOTIFF_NUM__3091__ProjCenterNorthingGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

GEOTIFF_NUM__3093__ProjScaleAtCenterGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

GEOTIFF_NUM__3094__ProjAzimuthAngleGeoKey=0.999600,45.309167,2546731.496000,-4354009.816000,337.255560,-86.000000

IMAGE__Z_RESOLUTION=0.000000

GEO__ModelTypeGeoKey=1

GEO__ProjectedCSTypeGeoKey=32767

GEO__PCSCitationGeoKey=IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection = Oblique Mercator (Hotine)|IMAGINE GeoTIFF Support

Copyright 1991 - 2001 by ERDAS, Inc. All Rights Reserved

@(#)$RCSfile: egtf.c $ $Revision: 1.7 $ $Date: 2001/09/27 14:59:29EDT $

Projection Name = Oblique Mercator (Hotine)

Units = meters

GeoTIFF Units = meters|

IMAGE__Z_ORIGIN=0.000000

IMAGE__COMPRESSION_VERSION=1,6,1

IMAGE__TARGET_COMPRESSION_RATIO=20.000000

IMAGE__COMPRESSION_NLEV=4

IMAGE__COMPRESSION_WEIGHT=4.000000

IMAGE__COMPRESSION_GAMMA=2.000000

IMAGE__COMPRESSION_BLOCK_SIZE=512

IMAGE__CREATION_DATE=Wed Mar 10 17:11:37 2004

VERSION=MG2

IMAGE__BITS_PER_SAMPLE=8

IMAGE__WIDTH=6415

IMAGE__HEIGHT=8061

IMAGE__COLOR_SCHEME=0

IMAGE__DATA_TYPE=0

IMAGE__XY_ORIGIN=663145.000000,240035.000000

IMAGE__X_RESOLUTION=1.000000

IMAGE__Y_RESOLUTION=1.000000

Image Structure Metadata:

INTERLEAVE=PIXEL

Corner Coordinates:

Upper Left ( 663144.500, 240035.500)

Lower Left ( 663144.500, 231974.500)

Upper Right ( 669559.500, 240035.500)

Lower Right ( 669559.500, 231974.500)

Center ( 666352.000, 236005.000)

Band 1 Block=1024x128 Type=Byte, ColorInterp=Red

Minimum=0.000, Maximum=255.000, Mean=97.677, StdDev=49.163

Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

Band 2 Block=1024x128 Type=Byte, ColorInterp=Green

Minimum=0.000, Maximum=255.000, Mean=97.555, StdDev=49.371

Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue

Minimum=0.000, Maximum=255.000, Mean=97.402, StdDev=49.024

Overviews: 3208x4031, 1604x2016, 802x1008, 401x504

–Adam

grass-user mailing list

grass-user@lists.osgeo.org

http://lists.osgeo.org/mailman/listinfo/grass-user

grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

William_Kyngesburye · February 11, 2011, 8:39pm

Hmm, something I've seen for a while (Mac only), and it didn't improve with the latest MrSID DSDK, is that there are errors in retrieving GCS info from a sid with GDAL - some parameters are reported wrong or dropped. It hasn't been a problem for me since I usually know the projection of the SID and can override it in the conversion.

It looks like the MrSID DSDK tool mrsiddecode gets it right.

On Feb 11, 2011, at 1:08 PM, Adam Dershowitz, Ph.D., P.E. wrote:

I am also working on OS X. But, I am using the MrSID plugin, from William Kyngesburye. So it seems that both gdalwarp and r.in.gdal are willing to import it. But, maybe I will try to search around for the apps that you are referring to, to see if they might do a better job on the conversion.
Are you suggesting that the MrSID plugin that I have is working, but is loosing reference data along the way?

--Adam

On Feb 11, 2011, at 10:47 AM, Stuart Edwards wrote:

MrSID files are a bit tricky. This is a proprietary format owned by LizardTech (see http://www.gdal.org/frmt_mrsid.html). That being said, they freely provide decode tools at their website (http://www.lizardtech.com/). Working on OS X, I use their Raster_DSDK in a little command line routine to convert the file to a geotiff (tifg) and then import it into GRASS with r.in.gdal. Note that the MrSID format is very efficient and boils a quad down to about 1.5 Gb. This will become 3 or 4 times greater as a tif so make sure you have plenty of room to store it and process it. I use an old version of Expressview (a LizardTech viewer that used to be available for OS X) to look at the image and determine the coordinates of the bit that I really want in order to minimize the resulting file sizes.

AFAIK there are no open source apps that deal with MrSID directly (for those that do, you have to install your own version of the decoder first)

Stu

-----
William Kyngesburye <kyngchaos*at*kyngchaos*dot*com>
http://www.kyngchaos.com/

First Pogril: Why is life like sticking your head in a bucket filled with hyena offal?
Second Pogril: I don't know. Why IS life like sticking your head in a bucket filled with hyena offal?
First Pogril: I don't know either. Wretched, isn't it?

-HitchHiker's Guide to the Galaxy

Adam_Dershowitz · February 11, 2011, 9:56pm

On Feb 11, 2011, at 12:14 PM, Stuart Edwards wrote:

I’m not familiar with the Kyngesburye plugin and haven’t used gdalwarp, but here:

I want to import it into a UTM zone 17 project. So I did the following:

gdalwarp -t_srs ‘+proj=utm +zone=17 +datum=WGS84’ input.sid output_utm.sid

it looks like you are just trying to create a new .sid file in zone 17 (output_utm.sid). From the meta file at the MI source, I think the file is already projected in utm - probably in zone 17, and with units of ‘meters’. So what you maybe intended was to just create a .tiffg from the .sid so that it can be imported into GRASS.

In mrsiddecode (the app that does all the work inside the plugin - I assume) I would just issue the command

./mrsiddecode -i Hancock.sid -o Hancock_south.tif -of tifg -ulxy 0 0 -lrxy 71145 29160

where

-i - input
-o - output
-of - output file type
ulxy - upper left x,y
lrxy - lower right x,y

and you get a nice geotiff that r.in.gdal will recognize - already in utm zone 17 most likely. Actually, in my case, the OH mrsids are projected in state plane coordinates but you get the picture.

To get the mrsiddecode app you must register at

http://developer.lizardtech.com/

so you can access to the Mac version of the decoder and mrsidinfo download page. And be slightly comfortable on the command line.

Stu

Nope, that’s not doing it. I downloaded mrsiddecode and ran it. It nicely generates a geotiff. But, if I try to import that I still get an error because the project is UTM, while the image is Oblique Mercator with omerc projection. I tried to use gdalwarp to go from the geotiff to utm, but I get an error:

gdalwarp -t_srs ‘+proj=utm +zone=17 +datum=WGS84’ oak_grove_sw.tifg oak_grove_sw_utm.tifg
ERROR 1: latitude or longitude exceeded limits
ERROR 1: Too many points (441 out of 441) failed to transform,
unable to compute output bounds.

If I run gdalinfo on the geotiff I get it looks to me like the projection has lat long near 0. So, I don’t think that the problem is a MrSID problem. But instead there is just something wrong with how I am doing the projection, or something is wrong with the file itself.

Here is the output of gdalinfo oak_grove_sw.tifg:

Driver: GTiff/GeoTIFF
Files: oak_grove_sw.tifg
Size is 6415, 8061
Coordinate System is:
PROJCS[“Projection = Oblique Mercator (Hotine)”,
GEOGCS[“NAD83”,
DATUM[“North_American_Datum_1983”,
SPHEROID[“GRS 1980”,6378137,298.2572221010002,
AUTHORITY[“EPSG”,“7019”]],
AUTHORITY[“EPSG”,“6269”]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433],
AUTHORITY[“EPSG”,“4269”]],
PROJECTION[“Hotine_Oblique_Mercator”],
PARAMETER[“latitude_of_center”,0.9996],
PARAMETER[“longitude_of_center”,0.9996],
PARAMETER[“azimuth”,0.9996],
PARAMETER[“rectified_grid_angle”,90],
PARAMETER[“scale_factor”,0.9996],
PARAMETER[“false_easting”,0.9996],
PARAMETER[“false_northing”,0.9996],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,“9001”]]]
Origin = (663144.500000000000000,240035.500000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 663144.500, 240035.500) ( 6d57’27.03"E, 3d 9’14.09"N)
Lower Left ( 663144.500, 231974.500) ( 6d57’25.53"E, 3d 4’52.96"N)
Upper Right ( 669559.500, 240035.500) ( 7d 0’53.76"E, 3d 9’12.88"N)
Lower Right ( 669559.500, 231974.500) ( 7d 0’52.24"E, 3d 4’51.77"N)
Center ( 666352.000, 236005.000) ( 6d59’9.64"E, 3d 7’2.93"N)
Band 1 Block=6415x1 Type=Byte, ColorInterp=Red
Band 2 Block=6415x1 Type=Byte, ColorInterp=Green
Band 3 Block=6415x1 Type=Byte, ColorInterp=Blue

Adam_Dershowitz · February 11, 2011, 10:40pm

On Feb 11, 2011, at 9:23 AM, Adam Dershowitz, Ph.D., P.E. wrote:

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.
I downloaded some aerial imagery from here:
http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211--,00.html

I want to import it into a UTM zone 17 project. So I did the following:
gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' input.sid output_utm.sid

I googled around some and it looks like this is actually an issue with the Michigan projection in Proj4.

I found this link:
http://www.osgeo.org/pipermail/mapserver-users/2006-July/018292.html

And this seems to give reasonable results, based on the above:

gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' -s_srs '+proj=omerc +lat_0=45.309166667 +lonc=-86.0 +alpha=337.255555556 +k=0.9996 +x_0=499839.8337 +y_0=528600.2398 +ellps=GRS80 +datum=NAD83 +units=m' oakgrove_sw.sid oakgrove_sw_utm.tifg

Then it imports into grass fine. So, it seems like something is broken with that projection.

--Adam

rshepard · February 11, 2011, 11:40pm

On Fri, 11 Feb 2011, Adam Dershowitz, Ph.D., P.E. wrote:

I googled around some and it looks like this is actually an issue with the
Michigan projection in Proj4.

I read that all of Michigan, and not just Detroit, is in serious trouble
so it doesn't suprise me that it's map projection is also broken.

Glad you got something to work.

Carpe weekend,

Rich

Stu1 · February 12, 2011, 2:38pm

Glad you found a fix. I'll have to remember this in case I have to venture north into MI.

It would be interesting to import the tifg into a new location in GRASS letting the image set the projection, and then reproject it from there into your project location with r.proj. But I guess since it's a an issue with PROJ that you'd probably run into the same problem.

Stu

On Feb 11, 2011, at 5:40 PM, Adam Dershowitz, Ph.D., P.E. wrote:

On Feb 11, 2011, at 9:23 AM, Adam Dershowitz, Ph.D., P.E. wrote:

I am trying to import some data in GRASS and have run into a problem. I hope someone can give me a little insight to help out.
I downloaded some aerial imagery from here:
http://www.michigan.gov/dnr/0,1607,7-153-10371_14546-30211--,00.html

I want to import it into a UTM zone 17 project. So I did the following:
gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' input.sid output_utm.sid

I googled around some and it looks like this is actually an issue with the Michigan projection in Proj4.

I found this link:
http://www.osgeo.org/pipermail/mapserver-users/2006-July/018292.html

And this seems to give reasonable results, based on the above:

gdalwarp -t_srs '+proj=utm +zone=17 +datum=WGS84' -s_srs '+proj=omerc +lat_0=45.309166667 +lonc=-86.0 +alpha=337.255555556 +k=0.9996 +x_0=499839.8337 +y_0=528600.2398 +ellps=GRS80 +datum=NAD83 +units=m' oakgrove_sw.sid oakgrove_sw_utm.tifg

Then it imports into grass fine. So, it seems like something is broken with that projection.

--Adam_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user