Conversion from lat/long to PCS to raster

classic Classic list List threaded Threaded
12 messages Options
Reply | Threaded
Open this post in threaded view
|

Conversion from lat/long to PCS to raster

Fernando Cacciola
Hi,

I plan to use libgeotiff (never used it beefore) to convert from
lat/long (GCS) to raster pixels.

IIIUC I need to convert lat/long to northings/eastings first (GCS ->
PCS), then PCS to Raster space.

I can see in libgeotiff the function GTIFPCSToImage, for PCS->Raster,
and I can also see a GTIFProj4FromLatLong(), which appears to be for
GCS to PCS, but I'm not sure how to use it.

Or rather, I can't get it to work: with a test geotiff I have,
lat/long and northings/eastings end up the same when I call
GTIFProj4FromLatLong (which of course won't give the correct raster
coordinates since this geotiff is ModelTypeProjected)

I called GTIFGetProj4Defn(&defn) FWIW.

Is this supported and I'm missing something? Or I cannot do this with
libgeotiff alone (and need GDAL instead)?

TIA

--
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Cassanova, Bill
If the image you have already contains the geotiff tags it is much easier to use gdal directly.  If the image does not yet contain the geotiff tags you can add them using gdal_translate and then use gdal directly.

Below is an example of reading a geotiff image and extracting the lat-lon.

General instructions:
GDALDatasetH my_data_set = GDALOpen(path.c_str(), GA_ReadOnly);.

OGRSpatialReferenceH my_spatial_reference = OSRNewSpatialReference(NULL);

OSRSetWellKnownGeogCS(my_spatial_reference, "WGS84");

char *pszDstWKT = NULL;

OSRExportToWkt(my_spatial_reference, &pszDstWKT);

OSRDestroySpatialReference(my_spatial_reference);

void *my_transformer = GDALCreateGenImgProjTransformer(my_data_set,
                                                          NULL,
                                                          NULL,
                                                          pszDstWKT,
                                                          FALSE, 0.0, 0);
float original_lat = some_lat_value;
float original_lon = some_lon_value;

float lat = original_lat;
float lon = original_lon;

int result =  GDALGenImgProjTransform(my_transformer, TRUE, 1, &lon, &lat,
                              NULL, &panSuccess);

If( result == 1 )
{
  Std::cout << original_lon << " is " << lon << " in pixel space" << std::endl;
  Std::cout << original_lat << " is " << lat << " in pixel space" << std::endl;
}

Hope this helps,
Bill


-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Fernando Cacciola
Sent: Friday, January 21, 2011 6:58 AM
To: [hidden email]
Subject: [Geotiff] Conversion from lat/long to PCS to raster

Hi,

I plan to use libgeotiff (never used it beefore) to convert from
lat/long (GCS) to raster pixels.

IIIUC I need to convert lat/long to northings/eastings first (GCS ->
PCS), then PCS to Raster space.

I can see in libgeotiff the function GTIFPCSToImage, for PCS->Raster,
and I can also see a GTIFProj4FromLatLong(), which appears to be for
GCS to PCS, but I'm not sure how to use it.

Or rather, I can't get it to work: with a test geotiff I have,
lat/long and northings/eastings end up the same when I call
GTIFProj4FromLatLong (which of course won't give the correct raster
coordinates since this geotiff is ModelTypeProjected)

I called GTIFGetProj4Defn(&defn) FWIW.

Is this supported and I'm missing something? Or I cannot do this with
libgeotiff alone (and need GDAL instead)?

TIA

--
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Frank Warmerdam
In reply to this post by Fernando Cacciola
On 11-01-21 06:57 AM, Fernando Cacciola wrote:

> Hi,
>
> I plan to use libgeotiff (never used it beefore) to convert from
> lat/long (GCS) to raster pixels.
>
> IIIUC I need to convert lat/long to northings/eastings first (GCS ->
> PCS), then PCS to Raster space.
>
> I can see in libgeotiff the function GTIFPCSToImage, for PCS->Raster,
> and I can also see a GTIFProj4FromLatLong(), which appears to be for
> GCS to PCS, but I'm not sure how to use it.
>
> Or rather, I can't get it to work: with a test geotiff I have,
> lat/long and northings/eastings end up the same when I call
> GTIFProj4FromLatLong (which of course won't give the correct raster
> coordinates since this geotiff is ModelTypeProjected)
>
> I called GTIFGetProj4Defn(&defn) FWIW.
>
> Is this supported and I'm missing something? Or I cannot do this with
> libgeotiff alone (and need GDAL instead)?

Fernando,

Are you linked against the PROJ.4 library? Without out that the
GTIFFProj4FromLatLong/ToLatLong() functions are NOPs though they
should return FALSE indicating failure.

I would suggest looking at bin/listgeo.c for an example of using
GTIFFProj4ToLatLong().

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, [hidden email]
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | Geospatial Programmer for Rent

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Fernando Cacciola
In reply to this post by Cassanova, Bill
Hi Bill,

This looks exactly what I need! Thank you.

Quick question:

Looking at the tags I can see that this particular TIFF has "4269/NAD83" as GCS,
rather than WGS84.

I guess I can just pass that string to OSRSetWellKnownGeogCS()?

Or perhaps I should create the new spatial reference from the tags in the image?


Best

--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Fernando Cacciola
In reply to this post by Frank Warmerdam
Hi Frank,
On 1/21/2011 11:53 AM, Frank Warmerdam wrote:

>
> Fernando,
>
> Are you linked against the PROJ.4 library? Without out that the
> GTIFFProj4FromLatLong/ToLatLong() functions are NOPs though they
> should return FALSE indicating failure.
>
Ah... indeed I didn't. Small detail ;)

> I would suggest looking at bin/listgeo.c for an example of using
> GTIFFProj4ToLatLong().
>
OK.

I'll try this before I jumpo to quickly to using GDAL instead as Bill suggested.

Best
--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Cassanova, Bill
In reply to this post by Fernando Cacciola
Hi Fernando,

I think you have to pass the string "NAD83" or "EPSG:4269" instead.

I am going off this info I found on page:

http://www.gdal.org/ogr/classOGRSpatialReference.html#096b8dde4fd2eb475acd376060940b02


-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Fernando Cacciola
Sent: Friday, January 21, 2011 10:04 AM
To: [hidden email]
Subject: Re: [Geotiff] Conversion from lat/long to PCS to raster

Hi Bill,

This looks exactly what I need! Thank you.

Quick question:

Looking at the tags I can see that this particular TIFF has "4269/NAD83" as GCS,
rather than WGS84.

I guess I can just pass that string to OSRSetWellKnownGeogCS()?

Or perhaps I should create the new spatial reference from the tags in the image?


Best

--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

noelkhan
In reply to this post by Fernando Cacciola
Fernando,

So long as you can get the tie-point and scale info, you probably could use gdal or libtiff to convert from PCS to Raster space (you would have to convert LAT/LON to PCS yourself; but that's just algebra).

Conceptually:

1. Convert LAT/LON to UTM to get (Xm,Ym) [see: http://www.rcn.montana.edu/resources/tools/coordinates.aspx]

2. Get the geotiff's tie-point and scales
    A: ModelTiePointTag: (I,J,0,  Xt,Yt,0)
         ModelPixelScaleTag (Sx,Sy,0)
    B: ModelTransformationTag (d=> Xt, h=>Yt)
         Raster coordinates: I=>0, J=>0

3. Compute Dx,Dy between point (Xt,Yt) and (Xm,Ym)

4. Scale Dx,Dy using Sx,Sy and add raster offset (I,J) to find (I',J')
    I' = I + (Dx * Sx)
    J' = J + (Dy * St)

Regards,
Noel
(who is also new to geotiffs)
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Fernando Cacciola
In reply to this post by Frank Warmerdam
Hi Frank,

I rebuilt libgeotiff to support PROJ.4
Since I am on Windows I had to tweak the build system by hand (as is, it did not
work, for a number of details I'll post about when I'm done)
It eventually all build OK but now I have a problem:

GTIFFProj4FromLatLong() gives me NANs

I suspect the problem is that I am not arranging for PROJ.4 to read any data, as
I am not sure how to do that exactly.

That is, what do I need to put where?

Any help appreciated.

TIA

--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Fernando Cacciola
In reply to this post by noelkhan
On 1/21/2011 5:39 PM, noelkhan wrote:
>
> Fernando,
>
> So long as you can get the tie-point and scale info, you probably could use
> gdal or libtiff to convert from PCS to Raster space

Right

> (you would have to
> convert LAT/LON to PCS yourself;

Which is precisely the part I'm missing now :)


> Conceptually:
>
> 1. Convert LAT/LON to UTM to get (Xm,Ym) [see:
> http://www.rcn.montana.edu/resources/tools/coordinates.aspx]
>
Thanks.. I needed thsi bit of theory in fact,


The second part I knew already (though I prefer to let libgeotiff do that for me)

Best

--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Fernando Cacciola
In reply to this post by Cassanova, Bill
Hi Bill,

> Hi Fernando,
>
> I think you have to pass the string "NAD83" or "EPSG:4269" instead.
>
> I am going off this info I found on page:
>
> http://www.gdal.org/ogr/classOGRSpatialReference.html#096b8dde4fd2eb475acd376060940b02
>

OK

FWIW, I couldn't succeed in getting this out using libgeotiff, but I succeeded
with gdal based on the code you posted!! Thank you :)


--
---
Fernando Cacciola
SciSoft Consulting, Founder
http://www.scisoft-consulting.com

_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

Kevan
In reply to this post by Fernando Cacciola
Fernando Cacciola <fernando.cacciola@...> writes:

>
> On 1/21/2011 5:39 PM, noelkhan wrote:
> >
> > Fernando,
> >
> > So long as you can get the tie-point and scale info, you probably could use
> > gdal or libtiff to convert from PCS to Raster space
>
> Right
>
> > (you would have to
> > convert LAT/LON to PCS yourself;
>
> Which is precisely the part I'm missing now :)
>
> > Conceptually:
> >
> > 1. Convert LAT/LON to UTM to get (Xm,Ym) [see:
> > http://www.rcn.montana.edu/resources/tools/coordinates.aspx]
> >
> Thanks.. I needed thsi bit of theory in fact,
>
> The second part I knew already (though I prefer to let libgeotiff do that for
me)
>
> Best
>

Hi,

Im just starting to use geotiffs as well. I have been having a hard time
figuring out how to really even start. I have all the lat/longs and have already
converted them into UTM. Im not sure what you are supposed to do for steps 2 to
4 mentioned earlier. If you could explain or post snippets of code with
explanation that would be great.

Kevan



_______________________________________________
Geotiff mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/geotiff
Reply | Threaded
Open this post in threaded view
|

Re: Conversion from lat/long to PCS to raster

noelkhan
DISCLAIMER: As noted above, I'm a newbie and its been a while since I played around with geoTiffs.

Kevan,

[Step 2]

CONCEPT: The idea behind step 2 (above) is to relate some coordinate on the image to some coordinate in the real world. Imagine dividing your image into a grid (of I-columns and J-rows). Define the top-left corner (I,J) of your image as (0,0).

Suppose you had an image depicting parcels in your neighborhood. Now suppose you superimposed this image onto Google Maps, such that the outline of parcels from the image matched the satellite imagry. While you were busy matching up parcel/fence lines, you paid no attention to the top-left corner of the image, but that's actually what we need. So what real-world coordinate (Xt,Yt) corresponds to the top-left corner (0,0) of your image?

[Steps 3 and 4]

CONCEPT: But what if your image doesn't quite match up with the satellite photo? What if the image looks like it needs to be stretched to properly overlay onto the satellite imagery? In this case, you need to scale your image along the image's x or y dimensions. Let's keep track of that stretching in (Sx, Sy).

Now imagine you have a diamond in the center of your image. Imagine that diamond corresponded to a baseball field on satellite imagery. Using the process described above, I could relate the left-most point of the diamond (I,J) to the left-most point of the baseball field (Xt,Yt). If the diamond is now adaquately superimposed over the baseball field, we're done. If not, we either (1) scale the image along a dimension or (2) change the tie point (Xt,Yt).

For example, suppose you had to scale the image to make the right-most point of the diamond match the right-most point of the baseball field. Using your image's grid, identify the deltaX between the left and right -most points of the diamond; call this Dx_img. Using Google Maps, similarly identify the difference between X-ordinates; let's call this Dx_wrd. Now find Sx: Dx_wrd = Dx_img x Sx. Similarly, find Sy.

PRACTICE: You determine the relationship between raster-space (your image) and the model-space (representation of the real world) by: (1) getting that information from the provider of your imagery or (2) trail-and-error; superimposing the image until it fits (described here). (3) Most people simply use specialized "georeferencing" software that allows them to identify points on the image and corresponding points in the model and let software superimpose handle all the details.

RECOMMENDED READING: GeoTiff Spec (http://trac.osgeo.org/geotiff/)

[Re Snipets]

Unfortunately I won't be much help here. I got involved with geotiffs while writing a library in VB that read tif tags and geotif keys. Anyhow, GDAL appears easy enough to use:

From http://lists.maptools.org/pipermail/mapserver-west/2005-June/000131.html 
> For multiple points (e.g. one for each corner) you just add another -gcp
> option with the values.  So, you might do something like:
>
> gdal_translate -gcp 0 0 -180 90  -gcp 0 300 -180 -90  -gcp 600 0 180 -90  -gcp
> 600 300 180 90 input.tif output.tif
>
> Then run it through gdalwarp to transform the image.  Hope that helps!
>
> Tyler
Regards,
Noel
(Newbie)