Hot answers tagged

106

The simplest way to transform coordinates in Python is pyproj, i.e. the Python interface to PROJ.4 library. In fact: from pyproj import Proj, transform inProj = Proj(init='epsg:3857') outProj = Proj(init='epsg:4326') x1,y1 = -11705274.6374,4826473.6922 x2,y2 = transform(inProj,outProj,x1,y1) print x2,y2 returns -105.150271116 39.7278572773 EDIT based on ...


100

selimnairb's answer is close but you wont have the headers unless you've installed libgdal-dev: sudo apt-get install libgdal-dev with that done, export CPLUS_INCLUDE_PATH=/usr/include/gdal export C_INCLUDE_PATH=/usr/include/gdal and then pip install GDAL The compilation ran to completion and I have GDAL in my virtual env. Phew! (edit 2018) Note: To ...


62

I do not know ArcPy, but I work with shapefiles and raster in Python for years For processing shapefiles in Python, there are many modules like osgeo/ogr, Fiona, Pysal or Pyshp (shpUtils is one of them and not the most used), and others, see Pypi: GIS and examples on gis.stackexchange and many examples on the Web (not only in English). Most of them are much ...


61

You can use the gdal.Dataset or gdal.Band ReadRaster method. See the GDAL and OGR API tutorials and the example below. ReadRaster does not use/require numpy, the return value is raw binary data and needs to be unpacked using the standard python struct module. An example: from osgeo import gdal,ogr import struct src_filename = '/tmp/test.tif' shp_filename ...


54

That's strange, as if people suddenly discovered the power of Python (without ArcPy which is just one Python module among others), see for example the question Visualize shapefile in Python: geospatial processing in Python has a very long history, much older than Arcpy (or arcgisscripting) -> no "mimic" the capabilities of ArcPy here, as Paul says, most ...


51

While shapely doesn't natively understand coordinate systems, shapely.ops.transform() can do that along with pyproj. If pyproj.Proj can understand your both of your coordinate systems, then it can be made into a function that shapely can transform with. From the shapely docs: from functools import partial import pyproj from shapely.ops import transform ...


51

EDIT I edited the answer for managing particular situations (due to specific angle values) and for not displaying the dotted lines when a round angle is defined. I propose a solution by only recurring to rule-based symbology and labeling. Before starting, I want to underline that I will focus the attention on the explanation of the minimal things to do for ...


46

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 ...


45

Well-known binary is a good binary exchange format that can be exchanged with plenty of GIS software, including Shapely and GDAL/OGR. This is a tiny example of the workflow with osgeo.ogr: from osgeo import ogr from shapely.geometry import Polygon # Here's an example Shapely geometry poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)]) # Now convert it to a ...


43

This can be done in far fewer lines of code src = gdal.Open(path goes here) ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform() lrx = ulx + (src.RasterXSize * xres) lry = uly + (src.RasterYSize * yres) ulx, uly is the upper left corner, lrx, lry is the lower right corner The osr library (part of gdal) can be used to transform the points to any ...


43

You could try shapely. They describe spatial relationships and it work on windows The spatial data model is accompanied by a group of natural language relationships between geometric objects – contains, intersects, overlaps, touches, etc. – and a theoretical framework for understanding them using the 3x3 matrix of the mutual intersections of ...


43

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': {'...


42

if you have python-gdal bindings: import numpy as np from osgeo import gdal ds = gdal.Open("mypic.tif") myarray = np.array(ds.GetRasterBand(1).ReadAsArray()) And you're done: myarray.shape (2610,4583) myarray.size 11961630 myarray array([[ nan, nan, nan, ..., 0.38068664, 0.37952521, 0.14506227], [ nan, nan, ...


41

I've been using "in_memory" quite a bit recently. It can be very useful, as it has the potential to dramatically increase processing speeds for certain tasks, however if you are working with very large datasets, it might cause your program to crash. You can use "in_memory" to define process outputs... often, if I am performing a task on a feature class, ...


41

Updated answer for Python 3 You should not nowadays encounter any issue while installing on Windows using (as long as pip is version 8 or more according to official Shapely doc and you are using Python 3): pip install shapely To check pip version, do : python -c "import pip;print(pip.__version__)" If you need to upgrade pip, just do : pip install -U ...


38

NumPy: NumPy is the fundamental package for scientific computing with Python. It contains among other things: a powerful N-dimensional array object sophisticated (broadcasting) functions tools for integrating C/C++ and Fortran code useful linear algebra, Fourier transform, and random number capabilities Besides its obvious scientific uses, NumPy can also ...


38

When you load a GeoJSON file using the json library, you get a dict that contains an entry features, which contains the list of features. Each feature in turn consists of a dict, which, among other things, contains an entry geometry. The geometry is a dict containing the entries type and coordinates. So you can traverse your GeoJSON file like this: import ...


37

Most methods to spline sequences of numbers will spline polygons. The trick is to make the splines "close up" smoothly at the endpoints. To do this, "wrap" the vertices around the ends. Then spline the x- and y-coordinates separately. Here is a working example in R. It uses the default cubic spline procedure available in the basic statistics package. ...


36

By default the site you linked to uses the Spatial Reference System EPSG 3857 (WGS84 Web Mercator). I found this information here. You can either specify another Spatial Reference System by entering the desired EPSG into the form under Spatial Reference or you can convert the returned coordinates with Python. For instance you can use the GDAL Python ...


35

Τhe header files cannot be found for some reason. Maybe you you are operating inside a Virtual Enviroment or they are not where they should be for some reason. In any case you can specify the include dirs when installing gdal via pip. first download python's gdal : pip install --no-install GDAL in later versions of pip (>= 9.0.0) pip install --no-install ...


35

Grab a copy of the ogr2ogr Python port, which is distributed with the GDAL source code download or can be found here: http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py Once you import that into your code, you can use it like this: import ogr2ogr def main(): #note: main is expecting sys.argv, where the first argument is the script name ...


35

If you've got an GDAL/OGR dev environment (headers, libs), you could radically simplify your code by using Fiona. To read features from a shapefile, add new attributes, and write them out as GeoJSON is just a handful of lines: import fiona import json features = [] crs = None with fiona.collection("docs/data/test_uk.shp", "r") as source: for feat in ...


35

NEW UPDATE QGIS3 can now be installed from osgeo's homebrew tap, which saves mucking about with python symlinks... i.e. brew tap osgeo/osgeo4mac brew install qgis Then type qgis in the terminal to launch, then pin the dock menu icon to simplify launching in future. Note that the installation console logs will also provide an example for how to symlink ...


34

I would recommend becoming familiar with the Python GDAL/OGR API to work with both vector and raster data. The easiest way to start using GDAL/OGR is via a python distribution such as python(x,y), Anaconda, or OSGeo4W. Further details on using GDAL for your specific tasks: Get Shapefile Fields and Types Get Projection Additionally, I would recommend the ...


32

It's a bit buried, but there is a second parameter to GetAttrValue() which returns the value at that ordinal. So I can do: In [1]: import osgeo.osr as osr In [2]: srs = osr.SpatialReference() In [3]: srs.SetFromUserInput("EPSG:27700") Out[3]: 0 In [4]: print srs PROJCS["OSGB 1936 / British National Grid", GEOGCS["OSGB 1936", DATUM["OSGB_1936",...


32

I didn't find a one-tool solution, but you can first use raster calculator to turn all values below a certain threshold to zero and then use gdal_translate with -a_nodata 0 to turn the 0 into nodata.


31

Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif. From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it ...


31

There are many modules to read shapefiles in Python, older than ArcPy, look at the Python Package Index (PyPi): shapefiles. There are also many examples in GIS SE (search for [Python] Fiona, for example) All can read the geometry, the fields and the projections. The older is osgeo (GDAL/OGR), look at the Python GDAL/OGR Cookbook for example another ...


30

I've designed Fiona (an OGR wrapper) to make this kind of processing simple. from fiona import collection import logging log = logging.getLogger() # A few functions to shift coords. They call eachother semi-recursively. def shiftCoords_Point(coords, delta): # delta is a (delta_x, delta_y [, delta_y]) tuple return tuple(c + d for c, d in zip(coords,...


Only top voted, non community-wiki answers of a minimum length are eligible