Re: [gdal-dev] ogr2ogr failing on postgis?

2016-04-10 Thread Even Rouault
Le dimanche 10 avril 2016 20:39:25, Paolo Cavallini a écrit :
> Il 10/04/2016 18:01, Even Rouault ha scritto:
> > Yes, this is expected and a rather common problem with shapefiles.
> 
> Thanks a lot for the explanation. QGIS therefore creates an incorrect
> prj, but I did not noticed it because it loaded it with the correct
> EPSG. Ticket upstream opened:
> http://hub.qgis.org/issues/14655

No, QGIS did it fine (actually I suspect it generates .prj file through OGR 
itself). ESRI .prj file are supposed *not* to contain EPSG codes unfortunately. 
This is the ESRI WKT variant. Which makes it hard to deal with them on the 
reading side. The issue is more on the OGR side.

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] ogr2ogr failing on postgis?

2016-04-10 Thread Paolo Cavallini
Il 10/04/2016 18:01, Even Rouault ha scritto:

> Yes, this is expected and a rather common problem with shapefiles.

Thanks a lot for the explanation. QGIS therefore creates an incorrect
prj, but I did not noticed it because it loaded it with the correct
EPSG. Ticket upstream opened:
http://hub.qgis.org/issues/14655
All the best, and thanks.
-- 
Paolo Cavallini - www.faunalia.eu
QGIS & PostGIS courses: http://www.faunalia.eu/training.html
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] ogr2ogr failing on postgis?

2016-04-10 Thread Even Rouault
Hi Paolo,

> Hi all,
> I have an ogr2ogr command, from shp to pg, that fails because it
> requests write access to spatial_ref_sys. An apparently identical
> command on another file runs smoothly. Any hint?
> Details and data documented here:
> http://hub.qgis.org/issues/14650

Yes, this is expected and a rather common problem with shapefiles. As the 
shapefile .prj doesn't contain any explicit EPSG code, the PostGIS driver tries 
to find a match in spatial_ref_sys based on an exact match on the WKT, but as 
the OGC WKT built from the .prj is different from the one in spatial_ref_sys 
(one of the difference is the absence of the AUTHORITY["EPSG","3003"] node of 
course !), it must fallback to creating an ad-hoc entry, hence the need for 
write access.

Workaround: provide explicit SRS, in that case -a_srs EPSG:3003 (as I can see 
in the .qpj)

A more "fuzzy" indentification of SRS could help for that.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] ogr2ogr failing on postgis?

2016-04-10 Thread Dmitry Baryshnikov

Hi Paolo,

It seems to me that your shape spatial reference not present in 
spatial_ref_sys and GDAL try to import it into spatial_ref_sys.

Your shape file have such SRS:
PROJCS["Monte_Mario_Italy_zone_1",GEOGCS["GCS_Monte 
Mario",DATUM["D_Monte_Mario",SPHEROID["International_1924",6378388,297]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",150],PARAMETER["false_northing",0],UNIT["Meter",1]]


If I execute such SQL: select * from spatial_ref_sys where srtext like 
'%Monte_Mario%' I received:


4265;"EPSG";4265;"GEOGCS["Monte 
Mario",DATUM["Monte_Mario",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6265"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745 
(...)"
4806;"EPSG";4806;"GEOGCS["Monte Mario 
(Rome)",DATUM["Monte_Mario_Rome",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6806"]],PRIMEM["Rome",12.452333,AUTHORITY["EPSG","8906"] 
(...)"
3003;"EPSG";3003;"PROJCS["Monte Mario / Italy zone 1",GEOGCS["Monte 
Mario",DATUM["Monte_Mario",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6265"]],PRIMEM["Greenwich",0,AUTHORITY[" 
(...)"
3004;"EPSG";3004;"PROJCS["Monte Mario / Italy zone 2",GEOGCS["Monte 
Mario",DATUM["Monte_Mario",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6265"]],PRIMEM["Greenwich",0,AUTHORITY[" 
(...)"
26591;"EPSG";26591;"PROJCS["Monte Mario (Rome) / Italy zone 1 
(deprecated)",GEOGCS["Monte Mario 
(Rome)",DATUM["Monte_Mario_Rome",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6806"]], 
(...)"
26592;"EPSG";26592;"PROJCS["Monte Mario (Rome) / Italy zone 2 
(deprecated)",GEOGCS["Monte Mario 
(Rome)",DATUM["Monte_Mario_Rome",SPHEROID["International 
1924",6378388,297,AUTHORITY["EPSG","7022"]],TOWGS84[-104.1,-49.1,-9.9,0.971,-2.917,0.714,-11.68],AUTHORITY["EPSG","6806"]], 
(...)"


I think you need to reproject shape file to EPSG: 3003 or defined it for 
shape file.


Best regards,
Dmitry

10.04.2016 18:23, Paolo Cavallini пишет:

Hi all,
I have an ogr2ogr command, from shp to pg, that fails because it
requests write access to spatial_ref_sys. An apparently identical
command on another file runs smoothly. Any hint?
Details and data documented here:
http://hub.qgis.org/issues/14650
Thanks.


___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

Re: [gdal-dev] Map algebra

2016-04-10 Thread Ari Jolma

08.04.2016, 19:28, Gregory, Matthew kirjoitti:

I would want any map algebra library to handle three types of issues when 
dealing with raster data:

   1) Handling mask or nodata values (e.g. an analysis mask should be
  able to be set on any algebra operation that retains or combines
  nodata values from its component operands)


GDAL supports nodata values and mask bands to identify cells that do not 
contain valid data. In arithmetics I think it is straightforward that  x 
+ nodata = nodata. In logical operations it depends on the 
interpretation of nodata. If it means unknown, then we are dealing with 
three-valued logics where for example unknown OR true = true. Think 
about for example a raster which shows polluted areas as a result of two 
studies (true, false, unknown). If one study shows a site polluted and 
the other has no data, then the result should be that the site is polluted.


I haven't yet got to the point in my code where I look at nodata and 
mask bands.




   2) Handling different spatial extents (e.g. the ability to specify
  an output extent -- coming from a specific raster or from the
  union/intersection of rasters)


This may be desirable in some cases. For example if we have a huge 
raster and want to study some parts of it. In general I would, however, 
leave this for other already existing GDAL tools, especially 
gdal_translate, whose functionality is available in APIs as well.




   3) Handling different cell sizes (e.g. the ability to specify an
  analysis cell size and resample/reproject(??) all operands based
  on this parameter)


I think using gdal_translate in this case as a first step makes even 
more sense.


Alex asked about how I iterate over blocks when there are two or more 
bands involved. The basic idea is that all methods operate on one band, 
and other bands (in my case if there is a second band) are adjusted to 
that. The algorithm is the following. "band" is a wrapper (struct) for a 
GDALRasterBand object and block cache. A block cache is simply an array 
of pointers to data, which are obtained from the GDALRasterBand object 
with ReadBlock() and which may be written back with WriteBlock().


band1 = gma_band_initialize(b1);
band2 = gma_band_initialize(b2);

while (iterate) {

block_index i;
for (i = blocks in band1) {

add block i to cache in band1;
block1 = get block i from band1;

update cache in band1 to allow desired focal distance;
update cache in band2 to allow desired focal distance;

call method specific callback with (band1, block1, 
band2, retval, arg);

make a note if the callback indicates a need for iteration;

based on the return value from the callback

return with error,
continue with next block from band1, or
write block1 to band1 and continue with next block 
from band1;


}
}
if iteration, then in some cases some band level things need to 
be done;

}
empty cache in band1;
empty cache in band2;
}

The cache update function reads all needed blocks into the cache and 
discards those that are not needed. If there were no discards, the cache 
would eventually contain the whole band.


In the method callback values from the band2 are obtained by first 
computing the global cell index and then the respective block and cell 
index within the block.


This may not be the best possible algorithm but to me it is pretty 
understandable, which is a goal in itself.


The code is at

https://github.com/ajolma/gdal/tree/trunk/gdal/map_algebra

The curse of dimensionality comes in at the point, where the above 
function is called. The method callback is a template function and it 
needs the datatype of each band that is given to it as an argument. (I 
see now that I have made also the above function a template function, 
which may not be needed).


Best,

Ari

___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev