Re: [gdal-dev] Storing Gdal datasets in Postgis
On Fri, Feb 26, 2010 at 2:51 PM, Ricardo Cezar Bonfim Rodrigues rikardoce...@msn.com wrote: Hi, I'm using Gdal to load dted files and get max elevations of a region, but I have constraints concerning speed. Inspite of loading geotifs generated from dted has reduced a lot the processing time, its still not enough to the constraints I have. I'm loading geotiffs which represents dted L1 files, reading its points using rasterIO and verifying the max elevation of a rectangle area given ( lat, lon, width). I've noticed the more expensive is the data reading, so I'm thinking about storing the dted (or geotiff) files in a Postgis. Would it increase the search speed? Does anyone know the best way to store dted files in a Postgis? I thought about storing geotiffs (generated from deted using gdal), but I dont now exactly how it would work. I know there are spatial queries in Postgis and I think it would be very useful in my case. I look forwarding to receving any tip of how to do this. Thanks in Advance, Ricardo Rodrigues Brazil Hi Ricardo, If you're planning to manage raster data in PostGIS, you should know we're working in a new extension (WKT Raster). Check this: http://trac.osgeo.org/postgis/wiki/WKTRaster It's under development yet, just like the GDAL WKT Raster driver: http://trac.osgeo.org/gdal/wiki/frmts_wtkraster.html Most of the spatial queries involving rasters are under development, I say, but it may be useful for you in the future... Keep an eye. Best regards, Jorge ___ 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] Storing Gdal datasets in Postgis
Hi Jorge, I've implemented a customized gdal reader using the method readBlock instead of using RasterIO and now it seems to be much faster and will propably be adequate to my problem. But thanks by the tips, I'll keep it in mind if I need it. Ricardo Rodrigues Brazil From: jorgearev...@gis4free.org Date: Mon, 1 Mar 2010 11:03:04 +0100 Subject: Re: [gdal-dev] Storing Gdal datasets in Postgis To: rikardoce...@msn.com CC: gdal-dev@lists.osgeo.org On Fri, Feb 26, 2010 at 2:51 PM, Ricardo Cezar Bonfim Rodrigues rikardoce...@msn.com wrote: Hi, I'm using Gdal to load dted files and get max elevations of a region, but I have constraints concerning speed. Inspite of loading geotifs generated from dted has reduced a lot the processing time, its still not enough to the constraints I have. I'm loading geotiffs which represents dted L1 files, reading its points using rasterIO and verifying the max elevation of a rectangle area given ( lat, lon, width). I've noticed the more expensive is the data reading, so I'm thinking about storing the dted (or geotiff) files in a Postgis. Would it increase the search speed? Does anyone know the best way to store dted files in a Postgis? I thought about storing geotiffs (generated from deted using gdal), but I dont now exactly how it would work. I know there are spatial queries in Postgis and I think it would be very useful in my case. I look forwarding to receving any tip of how to do this. Thanks in Advance, Ricardo Rodrigues Brazil Hi Ricardo, If you're planning to manage raster data in PostGIS, you should know we're working in a new extension (WKT Raster). Check this: http://trac.osgeo.org/postgis/wiki/WKTRaster It's under development yet, just like the GDAL WKT Raster driver: http://trac.osgeo.org/gdal/wiki/frmts_wtkraster.html Most of the spatial queries involving rasters are under development, I say, but it may be useful for you in the future... Keep an eye. Best regards,Jorge ___ 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
[gdal-dev] (no subject)
I am simply trying to open a raster image, manipulate in NumPy and then spit
it back out in GDAL.
I am having trouble setting null data...
I tried: ds.GetRasterBand(1).SetNoDataValue( - )
but nothing seems to work.
Any help is appreciated!
Below is my code:
#! /usr/bin/env python
import os, sys
use_numeric = True
try:
from osgeo import ogr, gdal
from osgeo.gdalconst import *
import numpy
os.chdir(r'L:\users\gv\numPyGDAL\ospy_data4\ospy_data4')
use_numeric = False
except ImportError:
import ogr, gdal
from gdalconst import *
import Numeric
os.chdir(r'L:\users\gv\numPyGDAL\ospy_data4\ospy_data4')
# register all of the drivers
gdal.AllRegister()
# open the image
ds = gdal.Open('sequest', GA_ReadOnly)
if ds is None:
print 'Could not open image'
sys.exit(1)
band = ds.GetRasterBand(1))# 1-based index
# read data and add the value to the string
data = band.ReadAsArray()
print data
dims = data.shape
nx=dims[1]
ny=dims[0]
nxarray=numpy.zeros(nx)
dst_format = 'HFA'
dst_datatype = gdal.GDT_Int16
dst_options = ['COMPRESS=LZW']
dst_file='sequest_final.img'
dst_xsize = nx
dst_ysize = ny
dst_nbands = 1
geoTransform = ds.GetGeoTransform()
proj = ds.GetProjection()
driver = gdal.GetDriverByName( dst_format )
dst_ds = driver.Create(dst_file, dst_xsize, dst_ysize, dst_nbands,
dst_datatype, dst_options)
for j in range(ny):
nxarray = data[j,:]
nxarray.shape = (1,nx)
print nxarray
dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j)
dst_ds.SetGeoTransform(geoTransform)
dst_ds.SetProjection(proj)
dst_ds = None
_
Hotmail: Free, trusted and rich email service.
http://clk.atdmt.com/GBL/go/201469228/direct/01/___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] (no subject)
Discourse Maps wrote: I am simply trying to open a raster image, manipulate in NumPy and then spit it back out in GDAL. I am having trouble setting null data... I tried: ds.GetRasterBand(1).SetNoDataValue( - ) but nothing seems to work. ... Dear Discourse Maps, Presumably you should be setting the nodata value on the first band of the dst_ds, not the source dataset. But I don't see any SetNoDataValue() call at all in your script. I would suggest adding it right after the Create() call. I would also note that the HFA driver does not have a COMPRESS=LZW creation option, though you can use [ 'COMPRESSED=YES' ] to use the grid compression scheme supported by the format. dst_format = 'HFA' dst_datatype = gdal.GDT_Int16 dst_options = ['COMPRESS=LZW'] dst_file='sequest_final.img' dst_xsize = nx dst_ysize = ny dst_nbands = 1 geoTransform = ds.GetGeoTransform() proj = ds.GetProjection() driver = gdal.GetDriverByName( dst_format ) dst_ds = driver.Create(dst_file, dst_xsize, dst_ysize, dst_nbands, dst_datatype, dst_options) for j in range(ny): nxarray = data[j,:] nxarray.shape = (1,nx) print nxarray dst_ds.GetRasterBand(1).WriteArray(nxarray,0,j) dst_ds.SetGeoTransform(geoTransform) dst_ds.SetProjection(proj) dst_ds = None Your email client has reformatted the script so it is hard to see the indentation. But I would advise strongly against doing SetGeoTransform() and SetProjection() many times on the same dataset. 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
[gdal-dev] Python Buffer Features
I'm writing some code that involves buffering a feature in a point layer and then doing a spatial filter with another set of layers. I've got everything working just fine but I'm having a hard time determining how the spatial filter works. Based on the documentation, it looks like all features that intersect the buffer will be included. Is this true? Is there any way to include only those features which are completely within the buffer? On a related note, how do I know what the units of the spatial filter are? Does it simply use the units of the original shapefile? Is there any way to use different units? Thanks, Spencer Gardner University of Wisconsin-Madison Department of Urban and Regional Planning ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Python Buffer Features
Spencer Gardner wrote: I'm writing some code that involves buffering a feature in a point layer and then doing a spatial filter with another set of layers. I've got everything working just fine but I'm having a hard time determining how the spatial filter works. Based on the documentation, it looks like all features that intersect the buffer will be included. Is this true? Is there any way to include only those features which are completely within the buffer? Spencer, The actual behavior may depend on whether you have built with GEOS support or not. Without it, you are guaranteed to get back all features which intersect the spatial filter geometry but there may be extras. In fact, in this situation the comparison is strictly done as a bounding box to bounding box comparison. If GDAL is built with GEOS support, then a final real intersect test is done using GEOS. Currently there is no way to request the spatial filter to return only geometries that are entirely within the spatial filter region. So you would need to do an extra test on the features you get to determine if they fall entirely within the filter geometry. I believe the OGRGeometry::Within() test is appropriate for this. On a related note, how do I know what the units of the spatial filter are? Does it simply use the units of the original shapefile? Is there any way to use different units? The spatial filter is interpreted as being in the same coordinate system as the layer it is being used with. If you are working with a shapefile then this will be the coordinate system of the shapefile. There is no provision for comparing or filtering with different units than the underlying layer. 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
[gdal-dev] Problems with gdalwarp re-projecting from lat-long to Lambert Conformal
Hello, I am having problems re-projecting an ArcAscii Grid from LatLong/WGS84 to Lambert Conformal. The source grid was derived from a GRASS mapset. Using r.out.gdal, a GRASS raster was exported to AAIGrid format using: r.out.gdal input=MODIS_input_raster format=AAIGrid type=UInt16 output=MODIS_output_latlong.asc The resulting grid has the following projection parameters: Driver: AAIGrid/Arc/Info ASCII Grid Size is 5039, 3244 Coordinate System is: GEOGCS[GCS_WGS_1984, DATUM[WGS_1984, SPHEROID[wgs84,6378137,298.257223563]], PRIMEM[Greenwich,0], UNIT[Degree,0.017453292519943295]] Origin = (-130.004289,47.007502) Pixel Size = (0.00833565,-0.00833565) Corner Coordinates: Upper Left (-130.0042888, 47.0075018) (130d 0'15.44W, 47d 0'27.01N) Lower Left (-130.0042888, 19.9666571) (130d 0'15.44W, 19d57'59.97N) Upper Right ( -88.0009546, 47.0075018) ( 88d 0'3.44W, 47d 0'27.01N) Lower Right ( -88.0009546, 19.9666571) ( 88d 0'3.44W, 19d57'59.97N) Center (-109.0026217, 33.4870795) (109d 0'9.44W, 33d29'13.49N) Band 1 Block=5039x1 Type=Int16, ColorInterp=Undefined NoData Value=255 When I try to re-project the exported AAIGrid to another projection, I use: gdalwarp -s_srs EPSG:4326 -t_srs '+proj=lcc +lat__1=33n +lat_2=45n +lon_0=97w' MODIS_output_latlong.asc MODIS_output_lambert.asc While the resultant projection information looks correct, Driver: GTiff/GeoTIFF Size is 5576, 4283 Coordinate System is: PROJCS[unnamed, GEOGCS[WGS 84, DATUM[unknown, SPHEROID[unnamed,6378137,298.2572235629972]], PRIMEM[Greenwich,0], UNIT[degree,0.0174532925199433]], PROJECTION[Lambert_Conformal_Conic_2SP], PARAMETER[standard_parallel_1,33], PARAMETER[standard_parallel_2,45], PARAMETER[latitude_of_origin,0], PARAMETER[central_meridian,-97], PARAMETER[false_easting,0], PARAMETER[false_northing,0], UNIT[metre,1, AUTHORITY[EPSG,9001]]] Origin = (-3540115.789204,5964029.040523) Pixel Size = (811.62081646,-811.62081646) Metadata: AREA_OR_POINT=Area Corner Coordinates: Upper Left (-3540115.789, 5964029.041) (142d23'11.19W, 42d57'29.68N) Lower Left (-3540115.789, 2487857.084) (128d 1'12.16W, 14d47'2.22N) Upper Right ( 985481.883, 5964029.041) ( 83d18'11.06W, 50d22'40.86N) Lower Right ( 985481.883, 2487857.084) ( 88d 2'31.57W, 19d32'55.81N) Center (-1277316.953, 4225943.062) (110d59'51.34W, 34d30'30.05N) Band 1 Block=5576x1 Type=Int16, ColorInterp=Gray The resultant text file has only the following on the top line: II* and possibly with lots of binary characters at the bottom. Thanks much for any help... Bill ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Georeferencing Geostationnary (+proj=geos) projection
Hello, I'm trying to put georeference information to a raster HDF file which doesn't containt projection nor geo-transformation or any GCP. HDF is not georeferenced, but contain lot of useful metadata, so I manage to get the projection, which is Geostationnary projection ( +proj=geos). However, I don't understand how to compute GeoTransformation or GCP to geo-reference properly my image: For classic flat geo-referenced image, computation between pixel and lat/lon can be done by gdal_translate -a_ullr. But it's not possible with GEOS projection since the image is not flat ! What is the correct method to put GeoTransformation or GCP to GeoTIFF ? Below some metadata of the file that could help. Thanks, Stéphane - Projection GeoReference=+proj=geos +h=3.58011e+07 +lon_0=73.9617 +ellps=WGS84 +x_0=-5632000 +y_0=5632000 Nominal Altitude (km)=42179.2 Nominal Central Point Coordinates (degrees)=0.349082, 73.9617 Nominal X_resolution (km)=2 Nominal Y_resolution (km)=2 Nominal False Northing (km)=5632 Nominal False Easting (km)=-5632 Nominal Pixel Stepping Angle (nanoradian)=223457.048379966 Nominal Line Stepping Angle (nanoradian)=223457.048379966 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Transform not working
I'm new to GDAL and am usingthe python bindings. I want to transform a newly
created shapefile from WGS84 to a Lambert Conical Conformal so that I can
calculate the length of the contained polylines. The script exits with no
errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate the
length directly from GDAL? Thanks.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
#~ feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1= 33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0')
coordtrans = osr.CoordinateTransformation(source, target)
geom.Transform(coordtrans)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] Re: Transform not working
Not much of a response going on here. In order to change the shapefile
projection from geographic to projected do I need to reopen the file,
get each feature, and transform the points individually? If so what
calls do I need to invoke?
On Mar 1, 12:54 pm, Jamie Lahowetz deadpic...@gmail.com wrote:
I'm new to GDAL and am usingthe python bindings. I want to transform a newly
created shapefile from WGS84 to a Lambert Conical Conformal so that I can
calculate the length of the contained polylines. The script exits with no
errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate the
length directly from GDAL? Thanks.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
#~ feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1= 33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0')
coordtrans = osr.CoordinateTransformation(source, target)
geom.Transform(coordtrans)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
___
gdal-dev mailing list
gdal-...@lists.osgeo.orghttp://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] Re: Transform not working
Jamie, Try transforming the line geometry before adding it to the feature.
That should do it.
chris
On Mon, Mar 1, 2010 at 4:24 PM, deadpickle deadpic...@gmail.com wrote:
Not much of a response going on here. In order to change the shapefile
projection from geographic to projected do I need to reopen the file,
get each feature, and transform the points individually? If so what
calls do I need to invoke?
On Mar 1, 12:54 pm, Jamie Lahowetz deadpic...@gmail.com wrote:
I'm new to GDAL and am usingthe python bindings. I want to transform a
newly
created shapefile from WGS84 to a Lambert Conical Conformal so that I can
calculate the length of the contained polylines. The script exits with no
errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate
the
length directly from GDAL? Thanks.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
#~ feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1= 33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0')
coordtrans = osr.CoordinateTransformation(source, target)
geom.Transform(coordtrans)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
___
gdal-dev mailing list
gdal-...@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
___
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] Re: Transform not working
Thanks for the response. I'm sorry I'm really new to this and am not sure
where to transform the line geometry, I know it should be before
layer.CreateFeature(feature) but after I add the points to the geometry. Is
it before or after feature.SetGeometry(line)? And is the right way to call
the transform
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1=33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0 +datum=NAD83')
coordtrans = osr.CoordinateTransformation(source, target)
line.Transform(coordtrans)?
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ source.SetFromUserInput('+proj=lcc +lat_1=33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0 +datum=NAD83')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ feature = layer.GetFeature(0)
#~ geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1=33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0 +datum=NAD83')
coordtrans = osr.CoordinateTransformation(source, target)
line.Transform(coordtrans)
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
On Mon, Mar 1, 2010 at 6:49 PM, Chris Garrard
garrard.chris+...@gmail.comgarrard.chris%2b...@gmail.com
wrote:
Jamie, Try transforming the line geometry before adding it to the feature.
That should do it.
chris
On Mon, Mar 1, 2010 at 4:24 PM, deadpickle deadpic...@gmail.com wrote:
Not much of a response going on here. In order to change the shapefile
projection from geographic to projected do I need to reopen the file,
get each feature, and transform the points individually? If so what
calls do I need to invoke?
On Mar 1, 12:54 pm, Jamie Lahowetz deadpic...@gmail.com wrote:
I'm new to GDAL and am usingthe python bindings. I want to transform a
newly
created shapefile from WGS84 to a Lambert Conical Conformal so that I
can
calculate the length of the contained polylines. The script exits with
no
errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate
the
length directly from GDAL? Thanks.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
Re: [gdal-dev] GDAL for BNG to Lat/Long, WGS84
Please? Anyone? I'm sure the response will be really quick... -- View this message in context: http://n2.nabble.com/gdal-dev-GDAL-for-BNG-to-Lat-Long-WGS84-tp4652000p4659130.html Sent from the GDAL - Dev mailing list archive at Nabble.com. ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL for BNG to Lat/Long, WGS84
scottd777 wrote: Hi, and thank heaven for GDAL! I was so happy to find this tool, but I am having some trouble figuring out the exact command I need for an image warping application. If there is someone who doesn't mind a newbie type question, I could really use the help! Essentially, I have imagery which is referenced with British National Grid. I am tyring to convert / warp this so that it can be used within Google Earth. I have reference points for the upper left and lower right of the image in Easting and Northing, and transformed those so that I got lat / long coordinates which line up very well within Google Earth, however it appears that the projection is off. I am hoping to use gdalwarp to convert the TIF file (referenced as mentioned), but I had some questions: - For the world file that I create for gdalwarp, should I include the Easting / Northing, or does this need to be lat / long? Scott, If the file is in the british national grid then the origin and pixel size should be in british national grid meters. - Can I just go straight to gdalwarp, or do I need to do some other translation first? If you have a TIFF with a world file you should be able to just proceed directly with the warp. - The command I am trying to use looks like this: gdalwarp -s_srs proj=utm +zone=33v +datum=OSGB36 -t_srs proj=latlong +datum=WGS84 filein.tif out.tif Is this the correct command to get from a BNG source ortho image to what would overlay correctly in something like Google Earth? I'm not google earth guru, but I had assumed that it expected imagery in the google mercator projection, not in geographic lat/long as you have targetted. I would instead use a command like: gdalwarp filein.tif out.tif -s_srs EPSG:27700 \ -t_srs +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgri...@null +wktext +no_defs I don't believe british national grid is the same as UTM 33, and I don't believe that a v in the zone indicator will mean anything to PROJ.4. EPSG:27700 is the british national grid. I have provided the full definition of the google mercator projection though shorter forms like EPSG:900913 or EPSG:3857 might work depending on the version of GDAL and supporting files you have. scottd777 wrote: Please? Anyone? I'm sure the response will be really quick... Not such an easy answer, and I'm not sure it will be everything you need. The problem with the popularity of Google Earth is that it has introduced a flood of new people with little GIS background all asking the same questions. 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] Re: Transform not working
deadpickle wrote:
Not much of a response going on here. In order to change the shapefile
projection from geographic to projected do I need to reopen the file,
get each feature, and transform the points individually? If so what
calls do I need to invoke?
On Mar 1, 12:54 pm, Jamie Lahowetz deadpic...@gmail.com wrote:
I'm new to GDAL and am usingthe python bindings. I want to transform a newly
created shapefile from WGS84 to a Lambert Conical Conformal so that I can
calculate the length of the contained polylines. The script exits with no
errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate the
length directly from GDAL? Thanks.
Jamie,
My suggestion in the following script would be to create the file with
the target srs instead of WGS84, and to transform the geometry
before even assigning it to the feature that will be written.
A few notes:
* There is no way with OGR of changing the coordinate system of
an existing layer. The only point at which it can be set is
layer creation time.
* The CreateFeature() method on the layer creates the feature in
the datastore (ie. shapefile) but there is no meaningful
connection maintained between the ogr.Feature object and the
feature in the file after this. Alterations to the feature
or it's geometry will have no effect on the file without doing
another layer.SetFeature() call.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
#~ feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1= 33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0')
coordtrans = osr.CoordinateTransformation(source, target)
geom.Transform(coordtrans)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
___
gdal-dev mailing list
gdal-...@lists.osgeo.orghttp://lists.osgeo.org/mailman/listinfo/gdal-dev
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
--
---+--
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
[gdal-dev] Re: Transform not working
I figured it out... with all of your help. I had the coordtrans =
osr.CoordinateTransformation switched. I wanted to input lat/lon (WGS84) and
transform it to lcc (source) and I also had to set the source prj to lcc. I
think you all said this stuff already, just took me a little while. Thanks
for all the help.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
target = osr.SpatialReference()
source.SetFromUserInput('+proj=lcc +lat_1=33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0 +datum=NAD83')
target.SetFromUserInput('WGS84')
coordtrans = osr.CoordinateTransformation(target, source)
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print float(value[3]),float(value[2])
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ trans from lat/lon input to lcc
line.Transform(coordtrans)
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
On Mon, Mar 1, 2010 at 12:54 PM, Jamie Lahowetz deadpic...@gmail.comwrote:
I'm new to GDAL and am usingthe python bindings. I want to transform a
newly created shapefile from WGS84 to a Lambert Conical Conformal so that I
can calculate the length of the contained polylines. The script exits with
no errors but when I look at the shapefile in arcCatalog I see that the
projection is still WGS84. Any ideas? Also, is there a way to calculate the
length directly from GDAL? Thanks.
import os
try:
from osgeo import ogr
from osgeo import osr
except:
import ogr
import osr
#~ start driver
driver = ogr.GetDriverByName('ESRI Shapefile')
#~ set workspace
ws = 'C:\\Users\\deadpickle\\
Desktop\\case study\\verify\\'
os.chdir(ws)
#~ check if it exists
if os.path.exists('test.shp'):
driver.DeleteDataSource('test.shp')
#~ create a new shp
ds = driver.CreateDataSource('test.shp')
#~ create spatref for new shp
source = osr.SpatialReference()
source.SetFromUserInput('WGS84')
#~ create a new data layer
layer = ds.CreateLayer('test',geom_type=ogr.wkbLineString,srs=source)
#~ add an id field
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
layer.CreateField(fieldDefn)
#~ create a geometry
line = ogr.Geometry(ogr.wkbLineString)
#~ get the def of the layer, scheme
featureDefn = layer.GetLayerDefn()
# create a new feature
feature = ogr.Feature(featureDefn)
#~ Read 2006track file
file = open(2006track.csv,r)
file.readline()
lnlist = file.readlines()
for lines in lnlist:
value = lines.split(,)
if float(value[1]) == 2:
print value
line.AddPoint(float(value[3]),float(value[2]))
#~ set the field 'ID'
feature.SetField('id', float(value[1]))
#~ set the geometry for this shp
feature.SetGeometry(line)
# add the feature to the output layer
layer.CreateFeature(feature)
#~ feature = layer.GetFeature(0)
geom = feature.GetGeometryRef()
#~ change projection to lcc
target = osr.SpatialReference()
target.ImportFromProj4('+proj=lcc +lat_1= 33.00 +lat_2=45.00
+lat_0=39.00 +lon_0=-96.00 +x_0=0.0 +y_0=0.0')
coordtrans = osr.CoordinateTransformation(source, target)
geom.Transform(coordtrans)
# destroy the geometry and feature and close the data source
ds.Destroy()
line.Destroy()
feature.Destroy()
file.close()
--
Jamie Ryan Lahowetz
University of Nebraska - Lincoln
Graduate Student - Geosciences
402.304.0766
jlahow...@gmail.com
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] Problems with gdalwarp re-projecting from lat-long to Lambert Conformal
On Mon, Mar 1, 2010 at 7:12 PM, Bill Hudspeth w...@unm.edu wrote: ... When I try to re-project the exported AAIGrid to another projection, I use: gdalwarp -s_srs EPSG:4326 -t_srs '+proj=lcc +lat__1=33n +lat_2=45n +lon_0=97w' MODIS_output_latlong.asc MODIS_output_lambert.asc Note that without -of format] you generate a GeoTIFF file... While the resultant projection information looks correct, Driver: GTiff/GeoTIFF ^^^ The resultant text file has only the following on the top line: ... it is binary, not a text file. check gdalwarp --formats for the available formats. Markus ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
