headers
Maciej Sieczka | 5 May 2008 20:47
Picon
Favicon

Re: SRS stuff broken in GDAL and PROJ

As the attachments were too big for the ML I'm re-sending the message
with the list of control points and script attached only.

Frank Warmerdam pisze:
> Maciej Sieczka wrote:

>> Ouch. So SRS stuff is pretty much broken in GDAL and PROJ. I 
>> subscribed to the ticket and I'd be happy to help with testing.

> I disagree with the assessment that things are broken.  They are just
>  a little bent. Numerically the differences are insignificant.

Frank,

So I checked and it looks that, unfortunately, things really are broken.

I have compared the results of re-projection from ll/WGS84 to each EPSG
code common for PROJ 4.5 and 4.6, using proj4 strings as present in the
"epsg" file.

Test points were obtained from the EPSG dataset v6.15 [1] using the
following SQL query on the SRS's area-of-use boundaries:

SELECT coord_ref_sys_code, area_west_bound_lon, area_south_bound_lat
FROM epsg_coordinatereferencesystem, epsg_area WHERE
epsg_coordinatereferencesystem.area_of_use_code=epsg_area.area_code AND
epsg_coordinatereferencesystem.show_crs=1 AND
epsg_coordinatereferencesystem.deprecated=0;

This gave a list of points located in the SW corner of each SRS's area
(Continue reading)

Permalink | Reply | 0 comments |
headers
Frank Warmerdam | 5 May 2008 20:57
Picon
Favicon
Gravatar

Re: SRS stuff broken in GDAL and PROJ

Maciej Sieczka wrote:
> Out of the 2865 common EPSG codes, proj4 strings of 133 of them are
> different enough between the epsg files of PROJ 4.5 and 4.6 to result in
> re-projection errors of few centimeters most of the time, about a meter
> in several cases to as much as kilometers in few. I would call this a
> critical issue.

Maciej,

I am presuming that the ones with large differences are due to the fact
that PROJ 4.6.0 does *not* apply a change-of-ellipse shift when no
datum to wgs84 mapping information is available.

I would suggest you isolate 2-3 examples out of your list that you are
confident are significant issues and we can discuss them in more detail.
Please ensure that you have worked out whether the significant difference
is in the PROJ.4 string generated for the epsg file, or whether it is
the behavior of PROJ.4 between 4.5 and 4.6 that you think is responsible.

It is important to isolate particular issues and try to understand them
deeply.

As you likely know Andrey has committed a fix in the last day or two
for CPLAtof() and with this change we should be able to generate an
epsg init file without the slight round off differences.

Best regards,
--

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam <at> pobox.com
(Continue reading)

Permalink | Reply | 0 comments |
headers
Maciej Sieczka | 6 May 2008 08:00
Picon
Favicon

Re: SRS stuff broken in GDAL and PROJ

Frank Warmerdam pisze:
> Maciej Sieczka wrote:

>> Out of the 2865 common EPSG codes, proj4 strings of 133 of them are
>>  different enough between the epsg files of PROJ 4.5 and 4.6 to
>> result in re-projection errors of few centimeters most of the time,
>> about a meter in several cases to as much as kilometers in few. I
>> would call this a critical issue.

> I am presuming that the ones with large differences are due to the
> fact that PROJ 4.6.0 does *not* apply a change-of-ellipse shift when
> no datum to wgs84 mapping information is available.

Right. After making correction for that, 120 SRSs in the epsg file of
PROJ 4.6 still have errors significant enough to cause an error of few
decimeters, compared to epsg file of PROJ 4.5.

> I would suggest you isolate 2-3 examples out of your list that you
> are confident are significant issues and we can discuss them in more
> detail.

All 120 are significant. Please pick for testing the one you like. I'm
attaching a list of test points I used - 120pts.txt. They are in
ll/WGS84 coordinates. The file is organized as follows:

---
espg_no easting northing
...
---
(Continue reading)

Permalink | Reply | 0 comments |
headers
Aydin | 7 May 2008 08:23
Picon

To handle geographic info from HDF5 formated Satellite Data

I try to handle geographic coordinates (lat/lon) from a sattellite data in 
HDF5 format. You can find necessary information below.
Could you please tell me how can initialize proj4 program and how can obtain 
lat/lon values for one pixel.

detailed info is below:

geo_number_rows = 4500
geo_number_columns = 3165
geo_pixel_size_x = 0.560648
geo_pixel_size_y = -0.56215554
geo_row_offset = -4759.928
geo_column_offset = -1636.3636
geo_dim_pixel = KM, KM
geo_pixel_def = CENTRE
geo_product_corners = -25.0, 25.0, -25.0, 75.0, 45.0, 75.0, 45.0, 25.0
geo_product_center = 9.719999, 50.47374

projection_name = EQUIDISTANT_CYLINDRICAL
projection_proj4_params = +proj=eqc +lat_ts = 0 +lon_0=0.000000

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply | 6 comments |
headers
Frank Warmerdam | 7 May 2008 18:00
Picon
Favicon
Gravatar

Re: To handle geographic info from HDF5 formated Satellite Data

Aydin wrote:
> I try to handle geographic coordinates (lat/lon) from a sattellite data in 
> HDF5 format. You can find necessary information below.
> Could you please tell me how can initialize proj4 program and how can obtain 
> lat/lon values for one pixel.
> 
> detailed info is below:
> 
> geo_number_rows = 4500
> geo_number_columns = 3165
> geo_pixel_size_x = 0.560648
> geo_pixel_size_y = -0.56215554
> geo_row_offset = -4759.928
> geo_column_offset = -1636.3636
> geo_dim_pixel = KM, KM
> geo_pixel_def = CENTRE
> geo_product_corners = -25.0, 25.0, -25.0, 75.0, 45.0, 75.0, 45.0, 25.0
> geo_product_center = 9.719999, 50.47374
> 
> projection_name = EQUIDISTANT_CYLINDRICAL
> projection_proj4_params = +proj=eqc +lat_ts = 0 +lon_0=0.000000

Aydin,

PROJ.4 does not address transformation from pixel/line space to georeferenced
space.  Just projected space to geographic space (and back).  The info above
does not indicate an ellipsoid, but we might assume WGS84 as a guess.
If the corners are in long/lat ordering we can assume the top left corner
is 25W, 75N.  This can be converted to equidistant cylindrical coordinates
like this:
(Continue reading)

Permalink | Reply | 5 comments |
headers
Marco Pasetti | 7 May 2008 20:53
Picon

DCW Projection

Hi all,

I need to import DCW (http://www.maproom.psu.edu/dcw/dcw_about.shtml) data
in GRASS, and I need to create a location with the correct projection to
import them.
On the DCW web site I read as follows:

COORDINATE SYSTEM DESCRIPTION

The coordinate system is described according to ARC/INFO standards. For more
information, see the project command in the ARC/INFO Users Guide.

General Projection Information

Projection:	GEOGRAPHIC
		The coordinate pairs are measured by latitude/longitude and
have not been altered by projection algorithms.
		If you are using ARCVIEW, the projection MUST be GEOGRAPIC
in order to reproject.

Units:		DD
		Coordinate pairs are measured in decimal degrees.

Spheroid:	CLARKE1866
		The spheroid for the data is the Clark 1866.

I created a location with the following projection: +proj=longlat
+ellps=clrk66 +no_defs

Is it correct? What would you suggest me?
(Continue reading)

Permalink | Reply | 0 comments |
headers
Aydin | 8 May 2008 12:56
Picon

Re: To handle geographic info from HDF5 formated Satellite Data

Thank you for response.
if I use your suggestion, for example pixel (1000,1000);

longitude = -25 + 0.560648 * 1000 = 5581.48
latitude = 75 - 0.560648 * 1000 = -5531.48

So, it is not possible.

I do not have problem converting projected space to geographic space using 
proj lib. My problem is how can calculate the distance from the corner using 
below parameters than I can convert those to lat/lon.

Best Regards,
Aydin

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply | 0 comments |
headers
Aydın Gürol Ertürk | 8 May 2008 12:58
Picon

Re: To handle geographic info from HDF5 formated Satellite Data

Frank,

Thank you for response.
if I use your suggestion, for example pixel (1000,1000);

longitude = -25 + 0.560648 * 1000 = 5581.48
latitude = 75 - 0.560648 * 1000 = -5531.48

So, it is not possible.

I do not have problem converting projected space to geographic space using
proj lib. My problem is how can calculate the distance from the corner using
below parameters than I can convert those to lat/lon.

Best Regards,
Aydin



2008/5/7, Frank Warmerdam <warmerdam <at> pobox.com>:
Aydin wrote:
I try to handle geographic coordinates (lat/lon) from a sattellite data in HDF5 format. You can find necessary information below.
Could you please tell me how can initialize proj4 program and how can obtain lat/lon values for one pixel.

detailed info is below:

geo_number_rows = 4500
geo_number_columns = 3165
geo_pixel_size_x = 0.560648
geo_pixel_size_y = -0.56215554
geo_row_offset = -4759.928
geo_column_offset = -1636.3636
geo_dim_pixel = KM, KM
geo_pixel_def = CENTRE
geo_product_corners = -25.0, 25.0, -25.0, 75.0, 45.0, 75.0, 45.0, 25.0
geo_product_center = 9.719999, 50.47374

projection_name = EQUIDISTANT_CYLINDRICAL
projection_proj4_params = +proj=eqc +lat_ts = 0 +lon_0=0.000000

Aydin,

PROJ.4 does not address transformation from pixel/line space to georeferenced
space.  Just projected space to geographic space (and back).  The info above
does not indicate an ellipsoid, but we might assume WGS84 as a guess.
If the corners are in long/lat ordering we can assume the top left corner
is 25W, 75N.  This can be converted to equidistant cylindrical coordinates
like this:

proj +proj=eqc +lat_ts=0 +lon_0=0.000000 +ellps=WGS84
-25 75
-2782987.27     8348961.81

So, the top left corner is -2782987.27mE, 8348961.81mN.

Actually, treating this as a projection issue is rather silly since
the projection is eqc which is just a scaling of geographic space.
And since all the locations and pixel sizes are expressed in decimal
degrees you might as well just ignore that aspect, and treat this as
a WGS84 image with the given corners.

So, if the top left corner is -25,75 and the pixel size is
0.560648, 0.56215554 you could compute pixel (20in from left,15 down
from top) as:

longitude = -25 + 0.560548 * 20
latitude = 75 - 0.560648 * 15

etc

Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam <at> pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org


_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply | 0 comments |
headers
Aydin | 8 May 2008 13:21
Picon

Re: To handle geographic info from HDF5 formated Satellite Data

Frank Warmerdam <warmerdam <at> pobox.com> writes:
> 
> Aydin,
> 
> PROJ.4 does not address transformation from pixel/line space to georeferenced
> space.  Just projected space to geographic space (and back).  The info above
> does not indicate an ellipsoid, but we might assume WGS84 as a guess.
> If the corners are in long/lat ordering we can assume the top left corner
> is 25W, 75N.  This can be converted to equidistant cylindrical coordinates
> like this:
> 
> proj +proj=eqc +lat_ts=0 +lon_0=0.000000 +ellps=WGS84
> -25 75
> -2782987.27     8348961.81
> 
> So, the top left corner is -2782987.27mE, 8348961.81mN.
> 
> Actually, treating this as a projection issue is rather silly since
> the projection is eqc which is just a scaling of geographic space.
> And since all the locations and pixel sizes are expressed in decimal
> degrees you might as well just ignore that aspect, and treat this as
> a WGS84 image with the given corners.
> 
> So, if the top left corner is -25,75 and the pixel size is
> 0.560648, 0.56215554 you could compute pixel (20in from left,15 down
> from top) as:
> 
> longitude = -25 + 0.560548 * 20
> latitude = 75 - 0.560648 * 15
> 
> etc
> 
> Best regards,
Frank,

Thank you for response.
if I use your suggestion, for example pixel (1000,1000);

longitude = -25 + 0.560648 * 1000 = 5581.48
latitude = 75 - 0.560648 * 1000 = -5531.48

So, it is not possible.

I do not have problem converting projected space to geographic space using 
proj lib. My problem is how can calculate the distance from the corner using 
below parameters than I can convert those to lat/lon.

Best Regards,
Aydin

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply | 2 comments |
headers
Frank Warmerdam | 8 May 2008 14:56
Picon
Favicon
Gravatar

Re: Re: To handle geographic info from HDF5 formated Satellite Data

Aydin wrote:
> Frank,
> 
> Thank you for response.
> if I use your suggestion, for example pixel (1000,1000);
> 
> longitude = -25 + 0.560648 * 1000 = 5581.48
> latitude = 75 - 0.560648 * 1000 = -5531.48
> 
> So, it is not possible.
> 
> I do not have problem converting projected space to geographic space using 
> proj lib. My problem is how can calculate the distance from the corner using 
> below parameters than I can convert those to lat/lon.

Aydin,

I stand corrected - the pixel size turns out to be in km not decimal degrees.

So you will have to offset from the top left corner using:

x (easting) = -2782987.27 + 560.648 * pixel
y (northing) = 8348961.81 - 560.648 * line

Then reproject that point back to long/lat using

proj -I +proj=eqc +lat_ts=0 +lon_0=0.000000 +ellps=WGS84

Good luck,
--

-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam <at> pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org

_______________________________________________
Proj mailing list
Proj <at> lists.maptools.org
http://lists.maptools.org/mailman/listinfo/proj
Permalink | Reply | 1 comment |

Gmane