Thread Subject: help with ERDAS IMAGINE import

I am tying to move a project from GRASS GIS to MATLAB (out future modeling efforts will be done in MATLAB, therefore there is a lot of motivation to be able to produce our model input grids in MATLAB. Our input layers will be a combination of raster grids, shape files, and netCDF data

I am having trouble with the first raster that I am trying to get into MATLAB. I realize that MATLAB will not read ERDAS IMAGINE (.img) files. Therefore I exported the file from GRASS as both a .mat file and as a geotif. I will describe both attempts below.

First the geotiff. After reading through the docs I was able to use geotiffinfo and geotiffread to create the following in my workspace:

1. foresttype_info (1 x 1 structure that appears to describe all of the georeferencing information correctly)
2. foresttype (4130 x 5456 matrix of the raster data)
3. foresttype_R (1 x 1 spatialref.MapRasterReference)
4. bbox (2 x 2 matrix bounding box)
5. cmap (65535 x 3 colormap matrix)
6. refmat (3 x 2 matrix... not sure what this is)

However, when I try and display "foresttype" with >> geoshow(map_data, foresttype_R), I get an empty figure with x and y axes 0:1.

Of note however, I am able to load the geotif into map viewer and display the correct looking projection with what appears to be the original color map from GRASS GIS.

I have also tried simply loading the .mat file. This produces a matrix of the raster data named "map_data" and four 1 x 1 matrices named "map_eastern_edge, map_northern_edge, map_southern_edge, and map_western_edge" along with the filename called "map_name".

Any suggestions to help break the impasse here would be much appreciated. Thanks in advance.

> First the geotiff. After reading through the docs I was able to use geotiffinfo and geotiffread to create the following in my workspace:
>
> 1. foresttype_info (1 x 1 structure that appears to describe all of the georeferencing information correctly)
> 2. foresttype (4130 x 5456 matrix of the raster data)
> 3. foresttype_R (1 x 1 spatialref.MapRasterReference)
> 4. bbox (2 x 2 matrix bounding box)
> 5. cmap (65535 x 3 colormap matrix)
> 6. refmat (3 x 2 matrix... not sure what this is)

refmat stands for "referencing matrix." You can ignore this variable; referencing matrices have been superseded by the more complete and self-documenting map raster reference (and geographic raster reference) objects.

> However, when I try and display "foresttype" with >> geoshow(map_data, foresttype_R), I get an empty figure with x and y axes 0:1.

You need to use mapshow instead of geoshow, because your raster is referenced to an x-y system of projected map coordinates -- a projected coordinate system (PCS). This is evident because the class of foresttye_R is spatialref.MapRasterence. If the data were unprojected -- in a latitude-longitude Geographic Coordinate System (GCS), then you'd get back a spatialref.GeoRasterReference object and it would be appropriate to use geoshow. The names (based on map vs. geo) should help you remember this:

If you have a spatialref.MapRasterReference, use mapshow.

If you have a spatialref.GeoRasterReference, use geoshow (with the map axes and projection of your choice).

By the way, look at foresttype_info to learn the specifics of the projected coordinate system (PCS).

> Of note however, I am able to load the geotif into map viewer and display the correct looking projection with what appears to be the original color map from GRASS GIS.

This make sense, because mapview, like mapshow, is intended for data sets that are already in a projected, x-y system.

> I have also tried simply loading the .mat file. This produces a matrix of the raster data named "map_data" and four 1 x 1 matrices named "map_eastern_edge, map_northern_edge, map_southern_edge, and map_western_edge" along with the filename called "map_name".

I can't follow what's going on with this step.

Hth,

Rob Comer
Mapping Toolbox Development
MathWorks

Thanks for the reply Rob.

"Rob Comer" <rob.comer.nospam@mathworks.com> wrote in message
> refmat stands for "referencing matrix." You can ignore this variable; referencing matrices have been superseded by the more complete and self-documenting map raster reference (and geographic raster reference) objects.
>
> > However, when I try and display "foresttype" with >> geoshow(map_data, foresttype_R), I get an empty figure with x and y axes 0:1.
>
> You need to use mapshow instead of geoshow, because your raster is referenced to an x-y system of projected map coordinates -- a projected coordinate system (PCS). This is evident because the class of foresttye_R is spatialref.MapRasterence. If the data were unprojected -- in a latitude-longitude Geographic Coordinate System (GCS), then you'd get back a spatialref.GeoRasterReference object and it would be appropriate to use geoshow. The names (based on map vs. geo) should help you remember this:
>
> If you have a spatialref.MapRasterReference, use mapshow.

I seem to having trouble re-prjecting the map with mapshow(). The examples I am looking at in the mapping toolbox webinar and the online examples all seem to be using geoshow on lat lon files.


I thought the same basic approach of

1. creating a new figure
2. defining a new projection with axesm()
3. calling mapshow().

I seem to be missing some step however, because this approach does not yield a new projection in new figure window in the same way it does with geoshow and a lat lon file.

In addition, I have been able to overlay all of my rasters with multiple calls mapshow() while using "hold on", however, the shape file I read into matlab with shaperead(), does not align properly with the rasters. It appears to be placed approximately 100 km south and west of the rasters (btw in the same projection and does overlay the rasters correctly in the GIS where I exported them). Any suggestions


>
> If you have a spatialref.GeoRasterReference, use geoshow (with the map axes and projection of your choice).
>
> By the way, look at foresttype_info to learn the specifics of the projected coordinate system (PCS).
>
> > Of note however, I am able to load the geotif into map viewer and display the correct looking projection with what appears to be the original color map from GRASS GIS.
>
> This make sense, because mapview, like mapshow, is intended for data sets that are already in a projected, x-y system.
>

Understood. Thanks.

"Kirk" <kwythers.nospam@umn.edu> wrote in message > I seem to having trouble re-projecting the map with mapshow(). The examples I am looking at in the mapping toolbox webinar and the online examples all seem to be using geoshow on lat lon files.
>
>
> I thought the same basic approach of
>
> 1. creating a new figure
> 2. defining a new projection with axesm()
> 3. calling mapshow().
>
> I seem to be missing some step however, because this approach does not yield a new projection in new figure window in the same way it does with geoshow and a lat lon file.

The first 2 steps should indeed set up a map axes with a new projection in a new figure window. You should make sure that the projection parameters match the map coordinate system to which the raster is registered. If they do, then the raster should appear correctly positioned with respect to the latitude-longitude graticule (aka "latitude-longitude grid").

> In addition, I have been able to overlay all of my rasters with multiple calls mapshow() while using "hold on", however, the shape file I read into matlab with shaperead(), does not align properly with the rasters. It appears to be placed approximately 100 km south and west of the rasters (btw in the same projection and does overlay the rasters correctly in the GIS where I exported them). Any suggestions?

This should work as long as the X and Y coordinates in the shapefile are in the same coordinate system as the raster. It's impossible to tell what's wrong here without looking at the datasets. Is the same unit of length in use everywhere? Can you open a help request with MathWorks tech support?

-- Rob