43
See the OGR Projections tutorial and the OGRSpatialReference class. In particular, the GetAttrValue method.
Here's a worked example.
from osgeo import gdal,osr
ds=gdal.Open(r'SOMERASTER.TIF')
prj=ds.GetProjection()
print prj
srs=osr.SpatialReference(wkt=prj)
if srs.IsProjected:
print srs.GetAttrValue('projcs')
print srs.GetAttrValue('geogcs')
For my ...
37
I created a small (and relatively popular) module which goes off and hunts for ArcGIS on your PC. Once find it adds the correct paths to the environment so that you can import arcpy. The usage goes like this:
try:
import archook #The module which locates arcgis
archook.get_arcpy()
import arcpy
except ImportError:
# do whatever you do if ...
34
You're doing your installation wrong.
Instead of
pip install shapely
go to https://pypi.org/project/Shapely/#built-distributions to see you can download Windows wheels at https://www.lfd.uci.edu/~gohlke/pythonlibs/#shapely
And On the second link, click on the file Shapely‑1.6.4.post1‑cp35‑cp35m‑win_amd64.whl
Shapely‑1.6.4 is the version of Shapely, cp35‑...
33
1) read your shapefile with Fiona, PyShp, ogr or ...using the geo_interface protocol (GeoJSON):
with Fiona
import fiona
shape = fiona.open("my_shapefile.shp")
print shape.schema
{'geometry': 'LineString', 'properties': OrderedDict([(u'FID', 'float:11')])}
#first feature of the shapefile
first = shape.next()
print first # (GeoJSON format)
{'geometry': {'...
17
use QgsMapLayer::loadNamedStyle
uri = "/home/user/style.qml"
layer.loadNamedStyle(uri)
16
If an unhandled exception, such as an ImportError, occurs before the add-in classes are instantiated they will become unresponsive, be given a [Missing] label, and have a red symbol for their icon in the case of items on toolbars or in menus.
You can confirm whether an import error is happening by wrapping your import statement with an exception handler and ...
15
Instead of doing the reclassification as a double for loop described by dmh126, do it using np.where:
# reclassification
lista[np.where( lista < 200 )] = 1
lista[np.where((200 < lista) & (lista < 400)) ] = 2
lista[np.where((400 < lista) & (lista < 600)) ] = 3
lista[np.where((600 < lista) & (lista < 800)) ] = 4
lista[np....
15
Use rasterio of Sean Gillies. It can be easily combined with Fiona (read and write shapefiles) and shapely of the same author.
In the script rasterio_polygonize.py
the beginning is
import rasterio
from rasterio.features import shapes
mask = None
with rasterio.drivers():
with rasterio.open('a_raster') as src:
image = src.read(1) # first band
...
14
With PyQGIS
ras = QgsRasterLayer("raster.tif")
pixelSizeX= ras.rasterUnitsPerPixelX()
pixelSizeY = ras.rasterUnitsPerPixelY()
print pixelSizeX
2.11668210081
print pixelSizeY
2.11685012701
With GDAL
from osgeo import gdal
raster = gdal.Open('raster.tif')
gt =raster.GetGeoTransform()
print gt
(258012.37107330866, 2.11668210080698, 0.0, 163176.6385398821, 0....
13
You are on the right track and the geopandas GeoDataFrame is a good choice for rasterization over Fiona. Fiona is a great toolset, but I think that the DataFrame is better suited to shapefiles and geometries than nested dictionaries.
import geopandas as gpd
import rasterio
from rasterio import features
Set up your filenames
shp_fn = '...
12
Something like this should work:
def reverse(s):
items = s.split()
digs = ''.join(i for i in s if i.isdigit())
dr = digs[::-1]
return ' '.join(map(None, items)).replace(digs, dr)
>>> reverse('321 test')
'123 test'
@mnpeterson brought up a good point about assuming where the numbers are...My post above would string all digits ...
12
Ran into this problem myself. If you want the x and y as separate GeoDataFrame columns, then this works nicely:
gdf["x"] = gdf.centroid.map(lambda p: p.x)
gdf["y"] = gdf.centroid.map(lambda p: p.y)
Starting with GeoPandas 0.3.0, you can use the provided x and y properties instead:
gdf["x"] = gdf.centroid.x
gdf["y"] = gdf.centroid.y
10
First case: insert new polyline
You have to define the Polyline like so:
pLine = [QgsPoint(1,1), QgsPoint(2,1), QgsPoint(2,2)]
it is a list-of-QgsPoint. After you can put more QgsPoint into the list by append to insert the element at the end of the list:
pLine.append(QgsPoint(3,2))
or inserting some vertex point at a specific location of the list:
## ...
10
Python window is not an equivalent of the Python shell; hence, you won't be able to use the raw_input there. To implement the interactivity with the user, you may choose any of these alternatives:
build custom script tools with input parameters (via arcpy.GetParameterAsText());
build Python add-ins (which have text boxes to fill in);
use 3rd party Python ...
9
You should check what the scale and offset are for your file. This can be done as follows:
van_taken.header.scale
van_taken.header.offset
This almost looks like an overflow error to me. The lower case x, y, and z properties need to re-scale and re-offset the coordinates to store it as an integer (which is how LAS files store them). To be honest, setting ...
8
The API docs state
Memory data providerType (memory)
The memory data provider is used to construct in memory data, for
example scratch data or data generated from spatial operations such as
contouring. There is no inherent persistent storage of the data. The
data source uri is constructed. The url specifies the geometry type
("point", "...
8
This sort of question is better answered in StackOverflow but the answer is straight-forward enough so I'll give you a hint here.
Your date is not a date as far as Python is concerned but a division sum - which is the main reason why it doesn't work. Your code also won't give you the last four digits. You need '[-4:]' (yours gives everything except the ...
8
To iterate over the layers in the map, get a reference to the current map document and list the layers within it:
mxd = arcpy.mapping.MapDocument("CURRENT")
layers = arcpy.mapping.ListLayers(mxd):
for layer in layers:
# do something with 'layer'
To just get the number of selected features, try leveraging the Describe object's fidSet property:
desc = ...
8
You can actually deal with each part of a multipart polygon by creating a separate polygon object. Take a look at the following code.
import arcpy
from arcpy import env
env.workspace = "C:/Data/Exercise08"
fc = "Hawaii.shp"
for row in arcpy.da.SearchCursor(fc, ["OID@", "SHAPE@"]):
print("Feature {0}: ".format(row[0]))
partnum = 0
for part in ...
8
For one thing, you should generally avoid explicit looping like that in Python whenever possible. Psycopg2's cursor objects support the iterator protocol. This means you can iterate row by row over the results without needing to manually take care of indices.
Another thing is that you are calling the execute function many times inside that loop when it ...
8
In cursors, length is a read-only property.
I couldn't imagine what a predictable outcome would be for setting a new length of a line. Would it just extend the last point out in the bearing from the next to last point? What if it were multipart? Would it grow the entire polyline segment by segment?
8
Ok, I'm sorry to post a question and then answer it myself so quickly, but I found a nice set of course slides from Utah State University that has a lecture on opening raster image data with GDAL. For the record, here is the code I used to open the PRISM Climate Group datasets (which are in the EHdr format).
def ReadBilFile(bil):
import gdal
gdal....
8
Something like below should work.
vLayer = iface.activeLayer()
canvas = iface.mapCanvas()
extent = vLayer.extent()
canvas.setExtent(extent)
8
The following approach uses a Search Cursor and Python dictionary to perform the following workflow:
Select points within each polygon feature
Update dictionary with key (OID) and value (point count) for each iteration
Find max point value and corresponding OID and write to a text file
import arcpy, os
points = r'C:\temp\mytest\points.shp'
polys = r'C:\...
8
I can crush this down to 3 lines of code, no cursors required!
import arcpy
arcpy.SpatialJoin_analysis("Site", "points","in_memory/points_SpatialJoin", "JOIN_ONE_TO_MANY", "KEEP_ALL", "", "INTERSECT")
arcpy.Statistics_analysis("points_SpatialJoin", "in_memory/stats", "Join_Count SUM","Id")
Then simply sort the table to find the polygon with most points.
8
Leaving the rest below, but the main thing was accessing the geometry properly. If iterating over rows, e.g. for index, row in zones.iterrows(): you can simply use row.geometry.centroid.x and row.geometry.centroid.y. Geometry is a special column included in a GeoDataFrame, so every row has a geometry attribute.
You are accessing that attribute, which ...
8
If you want to use the Calculate Field tool (instead of an Update Cursor), what you're assigning to val needs to be an unevaluated python expression. That is, it needs to be identical to the string you would type into the Calculate Field tool if you were using the GUI version.
What you're currently assigning to val is "!FEDIRP! !FENAME! !FETYPE!".
Any of ...
8
1) With ogr, see Fast tip: Filtering features using OGR Python for example
from osgeo import ogr
from os import remove
in_file= "shapefileA.shp"
out_file = "shapefileB.shp"
in_ds = ogr.Open( in_file )
in_lyr = in_ds.GetLayerByIndex(0)
if exists(out_file):
remove(out_file)
driver_name = "ESRI Shapefile"
drv = ogr.GetDriverByName( driver_name )
out_ds = ...
8
Because you imported unicode_literals, you're passing a unicode literal to GDAL which is expecting a string literal.
So, explicitly cast the 'GTiff' arg to a str.
E.g.
>>> from __future__ import unicode_literals
>>> from osgeo import gdal
>>> gdal.GetDriverByName('GTiff')
Traceback (most recent call last):
File "<stdin>...
7
While you can do it with the standalone installed copy of QGIS I would highly recommend installing QGIS and QGIS-dev using the OSGeo4W installer. OSGeo is easier to update and lets you get the latest stuff.
That will install Python into C:\OSGeo4W\app\Python27 and let you run Python from the OSGeo4W shell. python.exe lives in C:\OSGeo4W\bin.
I have some ...
Only top voted, non community-wiki answers of a minimum length are eligible
