Re: [gdal-dev] Shapefile - Multipolygon with interior rings
On Mon, Aug 9, 2010 at 4:51 PM, Jay L. wrote: > When printing the point_count (in blue) variable I am able to see that my > multipolygon consists of two polygons, each with 2 rings (one interior and > one exterior). I iterate through these rings and attempt to get a point > count, which I use to iterate over each linearring. The point count for the > exterior ring runs without a problem and returns the proper number of > points. The point count on the interior ring returns an error. Jay, It isn't immediately obvious to me why you are seeing the problem you report, but you are using the variable geom_count for both the outer and inner loop and this will cause you severe problems in any case where a multipolygon contains a different number of polygons than it's polygon's have rings. Best regards, -- ---+-- I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush | Geospatial Programmer for Rent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] How does GDAL's NetCDF driver determine the spatial reference system
James,
I added polar stereographic support to the trunk in r18931. I believe I
added all the coordinate systems to be cf=1.4 compliant. I did not have a
dataset to test it. I did not backport these changes to 1.7 because I was a
new commiter. I will test your file tomorrow on my 1.8dev machine.
kss
#
Kyle Shannon
Physical Science Technician
RMRS Fire Sciences Lab
Fire, Fuels & Smoke - RWU 4405
5775 Highway 10 W.
Missoula, MT 59808
(406)829-6954
kshan...@fs.fed.us
#
On Mon, Aug 9, 2010 at 5:50 PM, James Hiebert wrote:
> Sure. The ticket can be found here:
> https://trac.osgeo.org/gdal/ticket/3715
>
> ~James
>
> On Mon, Aug 09, 2010 at 03:13:59PM -0400, Lucena, Ivan wrote:
> > James,
> >
> > Could do file a ticked and upload a sample file on
> http://trac.osgeo.org/gdal/ ?
> >
> > Regards,
> >
> > Ivan
> >
> > James Hiebert wrote:
> > > Hi All,
> > >
> > > I'm trying to use GDAL (version 1.7.2 on Gentoo built with hdf5 and
> netcdf
> > > support) to read NetCDF files, but I'm a little unclear on how the
> NetCDF
> > > driver attempts to determine the georeferencing.
> > >
> > > The driver page here:
> > > http://www.gdal.org/frmt_netcdf.html
> > > states that it will attempt to read the metadata "spatial_ref" for a
> Well
> > > Know Text definition of the spatial reference (I assume that it means
> that
> > > it will try to use a global attribute of the NetCDF file named
> spatial_ref).
> > >
> > > Since this would be convenient for my application, I attempted to use
> GDAL
> > > this way without success.
> > >
> > > + I added a global attribute named "spatial_ref" with the WKT
> definition of
> > > the spatial reference system.
> > > + gdalinfo on the file still reports "Coordinate System is `'"
> > > + opening the file with the python bindings and calling
> dst.GetProjection()
> > > returns a blank string (i.e. no projection information)
> > >
> > > The full listings of my spatial reference definition, the output of
> ncdump
> > > and gdalinfo and the output from the python session are listed below.
> Can
> > > anyone tell me whether I'm doing something wrong, or should I not count
> on
> > > georeferencing being supported by the NetCDF driver?
> > >
> > > ~James
> > >
> > > ja...@basalt ~ $ ncatted -a
> spatial_ref,global,a,c,'PROJCS["unnamed",GEOGCS["Normal Sphere
> > >
> (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]'
> > > narccap_ecpc-20c3m-tasmin-ncep.nc
> narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
> > >
> > >
> > > ja...@basalt ~ $ ncdump -h narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
> > > netcdf narccap_ecpc-20c3m-tasmin-ncep-with-srs {
> > > dimensions:
> > > xc = 123 ;
> > > yc = 104 ;
> > > time = UNLIMITED ; // (9496 currently)
> > > bnds = 2 ;
> > > variables:
> > > char polar_stereographic ;
> > > polar_stereographic:grid_mapping_name =
> > > "polar_stereographic" ;
> > > polar_stereographic:longitude_of_central_meridian = 263. ;
> > > polar_stereographic:straight_vertical_longitude_from_pole =
> > > 263. ;
> > > polar_stereographic:standard_parallel = 60. ;
> > > polar_stereographic:TrueGridLength_Latitude = 60. ;
> > > polar_stereographic:XGridLength_M = 5. ;
> > > polar_stereographic:YGridLength_M = 5. ;
> > > polar_stereographic:PoleOnPlane = "north" ;
> > > polar_stereographic:false_easting = 3925000. ;
> > > polar_stereographic:false_northing = 8325000. ;
> > > polar_stereographic:latitude_of_projection_origin = 90. ;
> > > double yc(yc) ;
> > > yc:long_name = "y-coordinate of Cartesian system" ;
> > > yc:standard_name = "projection_y_coordinate" ;
> > > yc:axis = "Y" ;
> > > yc:units = "m" ;
> > > double xc(xc) ;
> > > xc:long_name = "x-coordinate of Cartesian system" ;
> > > xc:standard_name = "projection_x_coordinate" ;
> > > xc:axis = "X" ;
> > > xc:units = "m" ;
> > > double lon(yc, xc) ;
> > > lon:units = "degrees_east" ;
> > > lon:long_name = "longitude" ;
> > > lon:standard_name = "longitude" ;
> > > lon:axis = "X" ;
> > > lon:NX_NumPntsOnParallel = 123 ;
> > > lon:NY_NumPntsOnMeridian = 104 ;
> > > lon:actual_range = 211.5412f, 316.2759f ;
> > > lon:XGridLength_M = 5. ;
> > > lon:YGridLength_M = 5. ;
> > > lon:PoleOnPlane = "north" ;
> > > double lat(yc, xc) ;
> > > lat:units = "degrees_north" ;
> > > lat:long_name = "latitu
Re: [gdal-dev] How does GDAL's NetCDF driver determine the spatial reference system
Sure. The ticket can be found here:
https://trac.osgeo.org/gdal/ticket/3715
~James
On Mon, Aug 09, 2010 at 03:13:59PM -0400, Lucena, Ivan wrote:
> James,
>
> Could do file a ticked and upload a sample file on
> http://trac.osgeo.org/gdal/ ?
>
> Regards,
>
> Ivan
>
> James Hiebert wrote:
> > Hi All,
> >
> > I'm trying to use GDAL (version 1.7.2 on Gentoo built with hdf5 and netcdf
> > support) to read NetCDF files, but I'm a little unclear on how the NetCDF
> > driver attempts to determine the georeferencing.
> >
> > The driver page here:
> > http://www.gdal.org/frmt_netcdf.html
> > states that it will attempt to read the metadata "spatial_ref" for a Well
> > Know Text definition of the spatial reference (I assume that it means that
> > it will try to use a global attribute of the NetCDF file named spatial_ref).
> >
> > Since this would be convenient for my application, I attempted to use GDAL
> > this way without success.
> >
> > + I added a global attribute named "spatial_ref" with the WKT definition of
> > the spatial reference system.
> > + gdalinfo on the file still reports "Coordinate System is `'"
> > + opening the file with the python bindings and calling dst.GetProjection()
> > returns a blank string (i.e. no projection information)
> >
> > The full listings of my spatial reference definition, the output of ncdump
> > and gdalinfo and the output from the python session are listed below. Can
> > anyone tell me whether I'm doing something wrong, or should I not count on
> > georeferencing being supported by the NetCDF driver?
> >
> > ~James
> >
> > ja...@basalt ~ $ ncatted -a
> > spatial_ref,global,a,c,'PROJCS["unnamed",GEOGCS["Normal Sphere
> > (r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]'
> > narccap_ecpc-20c3m-tasmin-ncep.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
> >
> >
> > ja...@basalt ~ $ ncdump -h narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
> > netcdf narccap_ecpc-20c3m-tasmin-ncep-with-srs {
> > dimensions:
> > xc = 123 ;
> > yc = 104 ;
> > time = UNLIMITED ; // (9496 currently)
> > bnds = 2 ;
> > variables:
> > char polar_stereographic ;
> > polar_stereographic:grid_mapping_name =
> > "polar_stereographic" ;
> > polar_stereographic:longitude_of_central_meridian = 263. ;
> > polar_stereographic:straight_vertical_longitude_from_pole =
> > 263. ;
> > polar_stereographic:standard_parallel = 60. ;
> > polar_stereographic:TrueGridLength_Latitude = 60. ;
> > polar_stereographic:XGridLength_M = 5. ;
> > polar_stereographic:YGridLength_M = 5. ;
> > polar_stereographic:PoleOnPlane = "north" ;
> > polar_stereographic:false_easting = 3925000. ;
> > polar_stereographic:false_northing = 8325000. ;
> > polar_stereographic:latitude_of_projection_origin = 90. ;
> > double yc(yc) ;
> > yc:long_name = "y-coordinate of Cartesian system" ;
> > yc:standard_name = "projection_y_coordinate" ;
> > yc:axis = "Y" ;
> > yc:units = "m" ;
> > double xc(xc) ;
> > xc:long_name = "x-coordinate of Cartesian system" ;
> > xc:standard_name = "projection_x_coordinate" ;
> > xc:axis = "X" ;
> > xc:units = "m" ;
> > double lon(yc, xc) ;
> > lon:units = "degrees_east" ;
> > lon:long_name = "longitude" ;
> > lon:standard_name = "longitude" ;
> > lon:axis = "X" ;
> > lon:NX_NumPntsOnParallel = 123 ;
> > lon:NY_NumPntsOnMeridian = 104 ;
> > lon:actual_range = 211.5412f, 316.2759f ;
> > lon:XGridLength_M = 5. ;
> > lon:YGridLength_M = 5. ;
> > lon:PoleOnPlane = "north" ;
> > double lat(yc, xc) ;
> > lat:units = "degrees_north" ;
> > lat:long_name = "latitude" ;
> > lat:standard_name = "latitude" ;
> > lat:axis = "Y" ;
> > lat:NX_NumPntsOnParallel = 123 ;
> > lat:NY_NumPntsOnMeridian = 104 ;
> > lat:actual_range = 21.23807f, 67.63758f ;
> > lat:XGridLength_M = 5. ;
> > lat:YGridLength_M = 5. ;
> > lat:PoleOnPlane = "north" ;
> > ...
> >
> > // global attributes:
> > ...
> > :history = "Mon Aug 9 13:44:50 2010: ncatted -a
> > spatial_ref,global,a,c,PROJCS[\"unnamed\",GEOGCS[\"Normal
> > Sphere
> >
> > (r=6370997)\",DATUM[\"unknown\",SPHEROID[\"sphere\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latit
[gdal-dev] Shapefile - Multipolygon with interior rings
Good afternoon, I am working on a program to convert a shapefile from ocentric to ographic latitude. The program reads an ESRI shapefile, creates a new shapefile, copying the attribute data, and performs a latitude "shift" for each point in the geometry. OGR is working wonderfully on all geometry types, except for multipolygons with donuts. Below is the code snippet. Above this the geometry type and feature attributes are read from the input shapefile. When printing the point_count (in blue) variable I am able to see that my multipolygon consists of two polygons, each with 2 rings (one interior and one exterior). I iterate through these rings and attempt to get a point count, which I use to iterate over each linearring. The point count for the exterior ring runs without a problem and returns the proper number of points. The point count on the interior ring returns an error. Stepping up one "level" I can print the geometry (green text) for the multipolygon and the description of 2 polygons (in proper WKT format with an exterior and interior ring) are written. It seems that once I am "into" the linearrings of a multipolygon OGR is only "seeing" the exterior ring. Thoughts? Thanks! Jay Code Snippet: elif geom_type == "MULTIPOLYGON": while feature_in is not None: #Create the output feature from the input geometry = feature_in.GetGeometryRef() #Dig into the input geometry geom_count = geometry.GetGeometryCount() #Pack the multipolygon lunch for later multipoly = ogr.Geometry(ogr.wkbMultiPolygon) #Iterate through the polygons for g in range(geom_count): geom = geometry.GetGeometryRef(g) geom_count = geom.GetGeometryCount() geom_name = geom.GetGeometryName() #For individual polygons in the multipolygon layer if geom_name == "LINEARRING": #Create a linear ring to populate with the points. ring = ogr.Geometry(ogr.wkbLinearRing) #Get a point count point_count = geom.GetPointCount() #Iterate through the points for h in range (point_count): pntx = geom.GetX(h) pnty = geom.GetY(h) pntz = geom.GetZ(h) #perform the conversion if pntz == 0.0: new_y = convert(pnty) ring.AddPoint_2D(pntx, new_y) elif pntz != 0.0: new_y = convert(pnty) ring.AddPoint(pntx, new_y, pntz) #Add the linearring poly.AddGeometry(ring) #Assign the spatial reference to the feature poly.AssignSpatialReference(feature_in_spatref) #For multipolygons in the multipolygon layer elif geom_name == "POLYGON": #Create a polygon to populate with linear ring(s) poly = ogr.Geometry(ogr.wkbPolygon) #Iterate through the polygons for h in range (geom_count): geom = geom.GetGeometryRef(h) #Check for multipolyons with donuts try: #Get a point count point_count = geom.GetPointCount() #Create a linear ring to populate with the points. ring = ogr.Geometry(ogr.wkbLinearRing) #Iterate through the points for j in range (point_count): pntx = geom.GetX(j) pnty = geom.GetY(j) pntz = geom.GetZ(j) #perform the conversion if pntz == 0.0: new_y = convert(pnty) ring.AddPoint_2D(pntx, new_y) elif pntz != 0.0: new_y = convert(pnty) ring.AddPoint(pntx, new_y, pntz) #Exit gracefully if multipolygons with donuts are encountered except: print "This shapefile is not supported." #Add the linearring poly.AddGeometry(ring) #Add the polygon to the multipolygon multipoly.AddGeometry(poly) #Set the spatial reference of the feature multipoly.AssignSpatialReference(feature_in_spatref) ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] How does GDAL's NetCDF driver determine the spatial reference system
James,
Could do file a ticked and upload a sample file on http://trac.osgeo.org/gdal/ ?
Regards,
Ivan
James Hiebert wrote:
Hi All,
I'm trying to use GDAL (version 1.7.2 on Gentoo built with hdf5 and netcdf
support) to read NetCDF files, but I'm a little unclear on how the NetCDF
driver attempts to determine the georeferencing.
The driver page here:
http://www.gdal.org/frmt_netcdf.html
states that it will attempt to read the metadata "spatial_ref" for a Well
Know Text definition of the spatial reference (I assume that it means that
it will try to use a global attribute of the NetCDF file named spatial_ref).
Since this would be convenient for my application, I attempted to use GDAL
this way without success.
+ I added a global attribute named "spatial_ref" with the WKT definition of
the spatial reference system.
+ gdalinfo on the file still reports "Coordinate System is `'"
+ opening the file with the python bindings and calling dst.GetProjection()
returns a blank string (i.e. no projection information)
The full listings of my spatial reference definition, the output of ncdump
and gdalinfo and the output from the python session are listed below. Can
anyone tell me whether I'm doing something wrong, or should I not count on
georeferencing being supported by the NetCDF driver?
~James
ja...@basalt ~ $ ncatted -a
spatial_ref,global,a,c,'PROJCS["unnamed",GEOGCS["Normal Sphere
(r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]'
narccap_ecpc-20c3m-tasmin-ncep.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
ja...@basalt ~ $ ncdump -h narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
netcdf narccap_ecpc-20c3m-tasmin-ncep-with-srs {
dimensions:
xc = 123 ;
yc = 104 ;
time = UNLIMITED ; // (9496 currently)
bnds = 2 ;
variables:
char polar_stereographic ;
polar_stereographic:grid_mapping_name =
"polar_stereographic" ;
polar_stereographic:longitude_of_central_meridian = 263. ;
polar_stereographic:straight_vertical_longitude_from_pole =
263. ;
polar_stereographic:standard_parallel = 60. ;
polar_stereographic:TrueGridLength_Latitude = 60. ;
polar_stereographic:XGridLength_M = 5. ;
polar_stereographic:YGridLength_M = 5. ;
polar_stereographic:PoleOnPlane = "north" ;
polar_stereographic:false_easting = 3925000. ;
polar_stereographic:false_northing = 8325000. ;
polar_stereographic:latitude_of_projection_origin = 90. ;
double yc(yc) ;
yc:long_name = "y-coordinate of Cartesian system" ;
yc:standard_name = "projection_y_coordinate" ;
yc:axis = "Y" ;
yc:units = "m" ;
double xc(xc) ;
xc:long_name = "x-coordinate of Cartesian system" ;
xc:standard_name = "projection_x_coordinate" ;
xc:axis = "X" ;
xc:units = "m" ;
double lon(yc, xc) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:axis = "X" ;
lon:NX_NumPntsOnParallel = 123 ;
lon:NY_NumPntsOnMeridian = 104 ;
lon:actual_range = 211.5412f, 316.2759f ;
lon:XGridLength_M = 5. ;
lon:YGridLength_M = 5. ;
lon:PoleOnPlane = "north" ;
double lat(yc, xc) ;
lat:units = "degrees_north" ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:axis = "Y" ;
lat:NX_NumPntsOnParallel = 123 ;
lat:NY_NumPntsOnMeridian = 104 ;
lat:actual_range = 21.23807f, 67.63758f ;
lat:XGridLength_M = 5. ;
lat:YGridLength_M = 5. ;
lat:PoleOnPlane = "north" ;
...
// global attributes:
...
:history = "Mon Aug 9 13:44:50 2010: ncatted -a
spatial_ref,global,a,c,PROJCS[\"unnamed\",GEOGCS[\"Normal
Sphere
(r=6370997)\",DATUM[\"unknown\",SPHEROID[\"sphere\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",263],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3925000],PARAMETER[\"false_northing\",8325000]]
narccap_ecpc-20c3m-tasmin-ncep.nc
narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc\n",
"Fri May 8 09:34:14 2009: ncrcat
tasmin_ECPC_19
[gdal-dev] How does GDAL's NetCDF driver determine the spatial reference system
Hi All,
I'm trying to use GDAL (version 1.7.2 on Gentoo built with hdf5 and netcdf
support) to read NetCDF files, but I'm a little unclear on how the NetCDF
driver attempts to determine the georeferencing.
The driver page here:
http://www.gdal.org/frmt_netcdf.html
states that it will attempt to read the metadata "spatial_ref" for a Well
Know Text definition of the spatial reference (I assume that it means that
it will try to use a global attribute of the NetCDF file named spatial_ref).
Since this would be convenient for my application, I attempted to use GDAL
this way without success.
+ I added a global attribute named "spatial_ref" with the WKT definition of
the spatial reference system.
+ gdalinfo on the file still reports "Coordinate System is `'"
+ opening the file with the python bindings and calling dst.GetProjection()
returns a blank string (i.e. no projection information)
The full listings of my spatial reference definition, the output of ncdump
and gdalinfo and the output from the python session are listed below. Can
anyone tell me whether I'm doing something wrong, or should I not count on
georeferencing being supported by the NetCDF driver?
~James
ja...@basalt ~ $ ncatted -a
spatial_ref,global,a,c,'PROJCS["unnamed",GEOGCS["Normal Sphere
(r=6370997)",DATUM["unknown",SPHEROID["sphere",6370997,0]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",60],PARAMETER["central_meridian",263],PARAMETER["scale_factor",1],PARAMETER["false_easting",3925000],PARAMETER["false_northing",8325000]]'
narccap_ecpc-20c3m-tasmin-ncep.nc narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
ja...@basalt ~ $ ncdump -h narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc
netcdf narccap_ecpc-20c3m-tasmin-ncep-with-srs {
dimensions:
xc = 123 ;
yc = 104 ;
time = UNLIMITED ; // (9496 currently)
bnds = 2 ;
variables:
char polar_stereographic ;
polar_stereographic:grid_mapping_name =
"polar_stereographic" ;
polar_stereographic:longitude_of_central_meridian = 263. ;
polar_stereographic:straight_vertical_longitude_from_pole =
263. ;
polar_stereographic:standard_parallel = 60. ;
polar_stereographic:TrueGridLength_Latitude = 60. ;
polar_stereographic:XGridLength_M = 5. ;
polar_stereographic:YGridLength_M = 5. ;
polar_stereographic:PoleOnPlane = "north" ;
polar_stereographic:false_easting = 3925000. ;
polar_stereographic:false_northing = 8325000. ;
polar_stereographic:latitude_of_projection_origin = 90. ;
double yc(yc) ;
yc:long_name = "y-coordinate of Cartesian system" ;
yc:standard_name = "projection_y_coordinate" ;
yc:axis = "Y" ;
yc:units = "m" ;
double xc(xc) ;
xc:long_name = "x-coordinate of Cartesian system" ;
xc:standard_name = "projection_x_coordinate" ;
xc:axis = "X" ;
xc:units = "m" ;
double lon(yc, xc) ;
lon:units = "degrees_east" ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
lon:axis = "X" ;
lon:NX_NumPntsOnParallel = 123 ;
lon:NY_NumPntsOnMeridian = 104 ;
lon:actual_range = 211.5412f, 316.2759f ;
lon:XGridLength_M = 5. ;
lon:YGridLength_M = 5. ;
lon:PoleOnPlane = "north" ;
double lat(yc, xc) ;
lat:units = "degrees_north" ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
lat:axis = "Y" ;
lat:NX_NumPntsOnParallel = 123 ;
lat:NY_NumPntsOnMeridian = 104 ;
lat:actual_range = 21.23807f, 67.63758f ;
lat:XGridLength_M = 5. ;
lat:YGridLength_M = 5. ;
lat:PoleOnPlane = "north" ;
...
// global attributes:
...
:history = "Mon Aug 9 13:44:50 2010: ncatted -a
spatial_ref,global,a,c,PROJCS[\"unnamed\",GEOGCS[\"Normal
Sphere
(r=6370997)\",DATUM[\"unknown\",SPHEROID[\"sphere\",6370997,0]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",263],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",3925000],PARAMETER[\"false_northing\",8325000]]
narccap_ecpc-20c3m-tasmin-ncep.nc
narccap_ecpc-20c3m-tasmin-ncep-with-srs.nc\n",
"Fri May 8 09:34:14 2009: ncrcat
tasmin_ECPC_1979010106.nc tasmin_ECPC_1981010106.nc
tasmin_ECPC_1986010106.nc tasmin_ECPC_1991010106.nc
tasmin_ECPC_1996010106.nc tasmin_ECPC_2
Re: [gdal-dev] Raster Strategies Examples
The steps you've listed below are correct for maximizing performance when viewing the raster at (or around) the native resolution of the data. However, as you've found out, the performance can still be poor when trying to view the entire dataset. The basic rule is that your application (mapserver?) should not have to open more than 4 to 8 raster files to render the view. With that in mind (and since disk space is not a problem), create "layers" of the raster dataset at different scales, say a factor of 2, until you have only one image. For example, if the resolution of your original data is 1 meter per pixel, your layers would look like the following. Note the first eight layers are handled by the internal overviews you created with gdaladdo. Pixel Map Scale: Size: 1:n - -- 1 m 4,000 2 8,000 4 16,000 8 32k 16 64 32 128 64 250 128 500 Additional External Layers 256 1m 512 2m 10244m 20488m 409616m To create the additional external layers, use GDAL's gdalbuiltvrt to create a VRT file of the original dataset. You can then use gdal_retile.py or Maptiler to build the layers. Best Regards, Brent Fraser Edi KARADUMI wrote: I have read many strategies for raster performance, but i still have problems with my case. The posts that i have read explain the strategie, but are not very detailed. Im new to mapserver so i have problems implementing them. My case is: -about 6000 tiles that form the map -170mb each tile -aproximately 1.2T of iamges -i have 10T HD and i think disk space is not a problem for me -tile format is TIF -tile size 9375x6250 pixels The strategy i have implemented and the problems i have -First i divided the tiles in 60 folders, to increase the performance in disk seek/read. -Than i transformed source files into internally tiles with the command gdal_translate -of GTiff -co "TILED=YES" original.tif tiled.tif -i used gdaladdo to add internal pyramids gdaladdo -r average 2 4 8 16 32 64 128 -Created a tileindex using gdaltindex -write_absolute_path MapAll.shp //server/Maps/Subfolder1/*.tif -Created a spatial index .qix file shptree MapAll.shp -than added the layer to the mapfile without the .shp extension so the application can use the .qix as you may know, i have very slow performance when i zoom out and im stucked here. As i have read i should make a copy of the tiles with reduced resolution. Merge the tiles together and use min/max scale to show different layers in different scales. the min/max scales i zoom in/out are 100/120. Now my questions are - How can i calculate the scale where i should create another layer of the tiles, or i shoud see it with some tests? - how much should be the resolution of the new layer? - is there any tools or program to merge the tiles? merging 6000 tiles with the gdalwarp by writing the command by myself is frustrating - how many tiles should i merge together to create the new layer? (how many tiles should have the new layer) i know that in each zoomscale its better to appear only one tile but i dont know how to calculate it - the tiles that i should merge are the originals or those with internal tiling and overviews? ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Motion: Adopt RFC 29: OGR Set Ignored Fields
On Mon, Aug 9, 2010 at 2:19 AM, Peter J Halls wrote: > Please correct me if I am wrong, but is not the OGR SQL engine correct at > present, ie conformant to the SQL standard, whilst this proposal would > render OGR SQL non-conformant? > > For example, > > SELECT * FROM tablename WHERE condition > > is the most common case and is the form currently delivered by OGR; > > however, > > SELECT column2,column3,column7,columnn FROM tablename WHERE condition > > must be the second most common. SQL provides for the specification of > columns to return (* / all or a list) but does not, so far as I am aware, > provide a construct by which one can specify columns to be omitted from the > results (which, I suppose, would equate to * EXCEPT FOR). > > Surely the goal is compliance with the standards? Peter, RDBMS based providers will construct the list of columns to select by scanining through the OGR field definitions, and including all those not marked as excluded in the select statement. The fact that the high level API function excludes fields rather than including them has no effect at the lower levels. It is just a way of ensuring that things are not excluded by accident. I see no problem. Best regards, -- ---+-- I set the clouds in motion - turn up | Frank Warmerdam, warmer...@pobox.com light and sound - activate the windows | http://pobox.com/~warmerdam and watch the world go round - Rush | Geospatial Programmer for Rent ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Motion: Adopt RFC 29: OGR Set Ignored Fields
On Sat, Aug 7, 2010 at 4:33 PM, Frank Warmerdam wrote:
>
> Re: the TestCapability(), I'd like to see this added too.
> I think this could be done without restarting the RFC
> process.
I've done minor edits to the RFC, stating:
The drivers supporting this method will return TRUE to OLCIgnoreFields
("IgnoreFields") capability.
@Ragi, Jason: I understand your concerns, but since there are pros and
cons for both "desired" and "ignored" fields, the discussion which one
is better can be endless. I'd like keep to keep the RFC as it is.
I'll proceed with updating the original patch to comply with the final
RFC specification and create a ticket in Trac soon.
Regards
Martin
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
