qgis - How do I get utm data converted into wgs84 all within Python27?

How do I get utm data converted into wgs84 all within Python27?

Does anyone know how to import the QGIS Lyon 'reproject' tool I to Python27?

I have a Python script that writes utm Into a text file but I require it to write wgs84 lats and longs instead. All I need (I think) is a module that outputs lats and longs before getting written to the file.

I tried to install Proj and matplotlib and other suggested solutions I found online but the imports seem to fail. I need guidance how to do this successfully but that's for another question.

I have QGIS Lyon installed and I am familiar with reprojecting layers into wgs84 using the reproject tool built in. I know that it works and therefore there is something on my pc that does what I want to achieve. I am hoping that I'm able to simply import the module that is used in QGIS reproject tool.

coordinate system - Dymaxion/Butterfly alternative projection for QGIS

I've been reading the Dymaxion and Butterfly projections are not open source, so can not be used with QGIS.

Are there alternatives to these available in QGIS?

I'm looking for a single projection so I can export relatively-accurate SVG shapes of countries around the world when I need them.

It's important these countries are relatively-close to the way they would appear on a globe.

qgis - how to get first and last point of MultiLineString-objects in qgis3

i would like to identify the start- and endpoints within a set of polylines and return the coordinates of both points finally.

so far i have been able to iterate the features and found out that ArcGIS's export created a MulitLineString with a number of vertices (where the number corresponds to the length of the line). using feature.geometry().vertexAt(0) (where feature is the current feature within a for loop) returns the first point of every feature and the print-result looks like that but i didn't find a way to set up a "for vertex in vertices"-style iteration to get the last point of every line. the python2.x-style from Can vector layer get start point and end point of line using PyQGIS with

geom = feature.geometry().asPolyline()

start_point = QgsPoint(geom[0])*  

end_point = QgsPoint(geom[-1])* 

doesn't work anymore.asPolyline() seems to return an empty list.

I've tried to iterate the vertices with vertices(self) → QgsVertexIterator but cant get any results. since i am pretty new to pyqgis the QGIS Python API seems really confusing to me.

Answer

If you're input is a multilinestring, then you'll need to handle that by either iterating through all the parts or (blindly!) just taking the first part alone.

For QGIS <= 3.4:

To take the first part:

multilinestring = feature.geometry().get()
first_part = multilinestring.geometryN(0)
# first_part will be a QgsLineString object

To get the first/last vertices:

first_vertex = first_part.pointN(0)
last_vertex = first_part.pointN(first_part.numPoints()-1)

For QGIS > 3.4:

The API is much easier here. You can use code like this and it'll work for BOTH linestring and multilinestring inputs:

for part in feature.geometry().get():
    first_vertex = part[0]
    last_vertex = part[-1]

How to add feature class from ESRI personal geodatabase (.mdb) in QGIS

How do you add a feature class in a ESRI personal gdb in qgis. I understand this can be done through the option add vector layer and select geodatabase. But what after that ? It doesn't allow me to browse through .

arcgis desktop - Storing two kinds of point geometry in shapefile or creating multipoint shapefiles?

Following How many interior and exterior rings can a polygon have in a standard ArcGIS shapefile?, can I have a mixture of two kinds of geometry in a shapefile of type point. Just like I have about shapefiles of type polygon or I should always create a shapefile of type multipoint in order to support multipoints in my layer?

My problem is as you see at How can I use ogrinfo to reach information about a .shp?, I can't use ogrinfo in order to understand that the geometry type of a multipoint shapefile will be known as wkbMultipoint or wkbPointwhen I use this line of code:

OGRwkbGeometryType GeometryType = poLayer ->GetGeomType();  

for a layer of type multipoint before programming.

• Will OGR know a layer of type multipoint as wkbpoint or wkbmultipoint?



• When we create a shapefile in ArcCatalog, we have point and multipoint, but not multipolyline or multipolygon for polyline and polygon.Why?

Answer

Quoting the shapefile specification (page 4):

All the non-Null shapes in a shapefile are required to be of the same shape type.

You may choose point OR multipoint OR polygon, but you cannot use multiple geometry types in a single shapefile. Multi-part geometries are permitted in all flavors of multipoint, polyline, and polygon shapefiles.

arcgis desktop - How is the default label field chosen by ArcMap?

I have 2 questions:

I know how to manually go to Layer Properties>Labels tab>Label Field and choose my label field. I also know that I could set up my primary label field and save the setting as a layer (*.lyr) file. This question, though, is more about how the label field is chosen by the program prior to user interaction.

First: When a new layer (shapefile, ArcSDE layer, personal/file geodatabase featureclass) is added to ArcMap, how is the default label field chosen? From my experience, it isn't chosen by field order, nor alphabetical by field. My best guess so far is that it might chose the first field of text/varchar (alphanumeric) type.

Second: Is there a way to define what the "primary label field" will be for a layer within the shapefile, sde, or gdb files, or are you stuck with saving layer files or letting the system pick for you.

Answer

According to this post, the default display field is chosen according to the following priorities:

1. First field of type Text whose name contains the word "name" (case-insensitive)


2. First field of type Text


3. First field of an integer type (Long or Short, presumably)


4. First field of any type

I don't think there is any way to specify the primary display field without using a layer or layer file, other than using the above logic to name/order your fields accordingly.

data - Locating shapefile of European railways for goods transportion?

We are trying to make a project about European goods trade so we need railway data in shapefile format.

More specifically, we are trying to find the potential of mine products being exported from Turkey to Europe by using Geomarketing tools.

Basically, we need a European railways layer for ArcGIS (excluding passenger-only railways).

remote sensing - What should come first-Radiometric calibration or co-registration?

For a pair of images,in the pre-processing of satellite images, which step should come first- Radiometric calibration or co-registration?

Answer

Suppose that we have two images that we want to co-register or one image that we want to register to earth:

First step is to remove the errors in each image both geometrically and radiometrically. Each image has some geometric errors due to:

• Earth rotation


• Scan time skew


• Aspect ratio


• Panoramic effect (bowtie error)


• Earth curvature

These errors will cause the pixels to drift during image acquisition and so will effect radiometric information. So when we are removing (transferring the pixels to their correct position in image) these geometric errors, we should do radiometric interpolation, too. Radiometric interpolation can be done through:

• Bilinear interpolation



• Bicubic interpolation


• Nearest Neighbor

Also if the two images have different sizes, we should resize them in this step through the above interpolation techniques.

Second step is to co-register (determining the mathematical transformation between two image) the images. This can be done through different ways. One of them is to register both images to earth. When both images are registered to earth (the same reference system), they'll be coregistered to each other.
Different mathematical models are used for registeration based on different factors including the type of the sensor that is used to acquire the images. One of them that is used in HRSI images is terrain independent RPC coefficients

Thus we always remove the errors in each image first (calibration of each image) and then co-register them. This is true for all kinds of images in remote sensing including PolSAR, InSAR, Hyperspectral and Multispectral images.

Two images are co-registered when both of them are free of errors

coordinate system - Waterman butterfly projection in Mapnik

Like the title says, how would someone configure Mapnik to use the Waterman butterfly projection ?

Otherwise, what other tools would be able to render using this projection ?

Answer

I don't think mapnik or proj4 are able to render that kind of projection.

According to that excellent post, Openlayers with protovis library would be able to render not exactly the Waterman projection but the Fuller projection (also called Dymaxion).

You even have an online example here.

pyqgis - Printing centered map from a QGIS project for each point in shapefile?

I need to produce on the order of 100 maps that are centered on each point of interest in a shapefile. I would like to prepare all the layers in a master QGIS project, and set up the composition for one point (so that printing 100 maps could be done manually, if need be).

I'd have something like the following layers:

• basemap


• points


• point buffer

And I would like to then automate printing to svg something like:

1. For each point in a shapefile


2. Center the map canvas on that point



3. Filter a buffer layer to only have that point


4. Print an svg with composition from the map composer

I'm reasonably certain I know how to do 1 & 2, but haven't found details on 3 & 4 on this site.

Answer

1. In the print composer, enable atlas and use point layer as atlas coverage layer.


2. Set the map item to be controlled by atlas, and choose the fixed scale


3. 
   

   Back in the QGIS's main window, for each layer that you want to filter according to a certain distance to the point use the rule based symbology and use the following rule

   
   
   

   within($geometry, buffer(@atlas_geometry, distance))

   
   

How to get feature for a given bbox from shapefile by ogr

I have some features point,line,polygon with shapefile format, now I want to get features for a given boundbox, is it possible to remove the features out of the bbox?

And for a polygon, I think it is necessary to close it.

I wonder if this is possible?

postgis - Finding points along a path

I am unsure of the best approach to this problem. I can think of some ways of doing it which I will list below but I am looking for the best and most efficient way of doing this.

• Given a variable number of points, create a "route" or line between the points.



• Create a bounding box (of 50km[arbitrary]) around the route.


• Return all points within this bounding box


• Points must be returned in order along the route (A->D)

Note: The radius does not have to be circles. A,B,C,D can use rectangular bounding boxes.

Data

• Table with points (x,y - technically lat/lon)


• Input points to create a route



• Distance from route to include points

We are using PostgreSQL (8.4) and PostGIS (1.3.6) and Python.

This is the solution I came up with. Thoughts? Ideas? Input?

• Create three polygon "tubes" (A->B, B->C, C->D)


• Filter points to only those in the polygons.


• Subtract B->C from A->B so the polygons do not overlap (no duplicates)


• Subtract C->D from B->C so the polygons do not overlap (no duplicates)



• For each polygon


• Calculate vector (A->B, B->C, C->D)


• Find furthest point opposite of vector in the polygon.


• Calculate distance from furthest point to each site in polygon.


• Return points in order of distance (smallest -> largest)


• Join all the points together (A->B points, B->C points, C->D points)

My approach was naive, after taking some time and looking at the PostGIS methods I came up with this single SQL call. Note its specific to our database but might help someone in the future. Im sure it could be improved as well.

    SELECT aerodrome_id

    FROM
      geo_aerodromecustom_view, 
      (SELECT ST_Buffer(
          ST_Transform(
              (SELECT ST_MakeLine(geom) 
              FROM 
                  (SELECT geom 
                  FROM geo_aerodromecustom_view
                  WHERE aerodrome_id IN ('CYOW','CYVR') 
                  ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END

                  ) sites
              ),
              900913
          ),
          50000) as geo
      ) as line
    WHERE ST_Intersects(ST_Transform(geom, 900913), line.geo)
    ORDER BY ST_Line_Locate_Point(
      ST_Transform(
          (SELECT ST_MakeLine(geom) 

           FROM 
              (SELECT geom 
               FROM geo_aerodromecustom_view
               WHERE aerodrome_id IN ('CYOW','CYVR') 
               ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END
              ) sites
          ),
          900913
      ),
      ST_Transform(geom, 900913)

    )

I had to project the data because I am using PostGIS 1.3.4 which doesnt support Geography type.

Basically what I am doing is I am using ST_MakeLine and a query to locate "aerodromes" and return their geometry.

I had to order them (using the CASE directive) so that the line would be connected in the right order.

I then project and buffer this line to create a Polygon that I can then use to see what other aerodromes intersect with the buffered polygon.

Using the unbuffered line (route) and the call ST_Line_Locate_Line I then order the discovered aerodromes as they appear appear along the path.

Answer

You need two postgis functions ST_Buffer and ST_Line_Locate_Point.

Access Violation in ArcObject Multi-threaded Application

I think I have figured it out. The most probable cause of exception in Multithreaded application may be:

When both threads trying to open same featureclass it breaks: coz one of the thread is already in the process of opening it and other tries to do the same, So it gives error like:

"Memory could not be read or write from protected memory"

This is my guess, but what wonders me is that why there is no separation of concern, even if I am opening a different workspace in each thread. May be Arcobjects internally, looks for same address space for any object in database. I have written few lines of code and tried to run it with different threads opening different object classes and this also breaks but

whenver, I deligate few thread to open same featureclass it breaks with the above errors.

Also, not releasing memory is the other cause for such errors as well.

Any suggestions?

Adding some code for more clarification:

private void ReadFeatureClass(string featureclassname)  
{    
    IWorkspaceFactory workspaceFactory = new SdeWorkspaceFactory();  

    IPropertySet connectionProperties = new PropertySet();  
    connectionProperties.SetProperty("SERVER",settings.Server); 

    connectionProperties.SetProperty("INSTANCE", Settings.Instance);  
    connectionProperties.SetProperty("USER", settings.User);  
    connectionProperties.SetProperty("PASSWORD", settings.Password);  
    connectionProperties.SetProperty("VERSION", settings.Version);  

    FeatureWorkspace featureWorkspace = workspaceFactory.Open(connectionProperties, 0) as IFeatureWorkspace;  

    ITable featureClass = featureWorkspace.OpenTable(featureClassName);  
}  

Above code breaks with the said exception most of the time.

@AndOne, I am using oracle spatial direct connect. Also, I created a new feature class and ran the same code against it, this also fails when threads are increased. There are few differences which I could figure out between these two featureclasses:

Any suggestions?

sql - PostGIS Intersection and Summarise Attributes

I am usually an ArcGIS desktop user, but keen to start using PostGIS more and I have a really big bit of processing to do. Not sure what functions to use, hopefully someone can help.

I have a polygon dataset (several million features), based on a type of landuse/ landcover classfication (20 categories). I have a number of regions in another dataset.

For each of the regions, I would like to know the area of each landcover classfication.

In ArcGIS (if it was a smaller dataset) I would imagine first adding the region to each of the polygons in the attribute table using a join. Then using "summarize" on the table by region and by landcover classification.

Not sure where to start doing this in PostGIS / SQL.

Update:

Wow thanks that has been a huge help.

It has been running a long time (44 hours!) and now I get:

NOTICE:  TopologyException: found non-noded intersection between LINESTRING (coords
edited) and LINESTRING (coords edited) at (position edited)
ERROR:  GEOS Intersection() threw an error!

********** Error **********

ERROR: GEOS Intersection() threw an error!
SQL state: XX000

I assume this is a problem in the original data - just a case of reviewing the original data or can I first check the topology some how for the whole data? Is there something about accepting certain errors / processing tolerances?

Answer

Assuming you have the following table layout

landcover(id,type,the_geom)
region(id,name,the_geom)

Area values per landcover type per region can be calculated using

SELECT r.id, r.name, r.the_geom, l.type, 
       Sum(ST_Area(ST_Intersection(l.the_geom,r.the_geom)))
FROM landcover l,
     region r

WHERE l.the_geom && r.the_geom
GROUP by r.id, r.name, r.the_geom, l.type

ST_Intersection is used to account for landcover polygons that are only partially within a region.

arcgis desktop - How to make a route event layer?

I'm using ArcGis 10.1 but I'm a beginner in this program.

I want to make a route event layer based on a route I have already defined. My question is in the third line when the program ask me for an input event table. I read a couple of explanations about this tool but any of them explains what table is this. Do I have to create an excel table? If yes, how does it have to look like?

Answer

You need to supply a table. This can be a .dbf file, a geodatabase table, or a sheet of an Excel file. I'd recommend exporting Excel/.dbf to a geodatabase table to make sure the field data types are converted properly.

After you've supplied the table, you will need to provide several fields (basically map your input table fields to the required in-fields). You can read what those fields mean here at the Make Route Event Layer (Linear Referencing) help page.

If you are new to Linear Referencing, consider going through ArcTutor tutorial which is shipped with your ArcGIS media. You can also download it from the Esri Customers Care portal.

qgis - Upright/Horizontal labels when labeling circular polygons at perimeter?

I have several polygons (circles) I want to label with their ID number. The label should be outside the circles so I use the positioning options "Use perimeter" and the tickbox "Next to the line (or similar, my qgis talk spanish...).

QGIS now automatically aligns the labels with the curvature of the line. Is there an option or a way to make them oriented horizontally?

pyqgis - How to configure QGIS (SAGA, GRASS algorithms) in Ubuntu

I have installed GRASS GIS following this install manual.

I have installed SAGA GIS following this install manual.

and finally I have installed QGIS 2.18.1 following the official details for install in Ubuntu Debian/Ubuntu.

The problem is I can't use algorithms from GRASS, SAGA and TAUDEM in the QGIS interface.

I have toolbox processing plugin and I have configure the options providers for all software (SAGA, GRASS) and not working.

Only GDAL tools are working in the plugins and QGIS algorithms. Standalone out of QGIS, SAGA and GRASS GIS work fine.

How to configure QGIS to work like work in Windows?

System details QGIS 2.18.1 Ubuntu 16.04.

GRASS ERRORS:

SAGA GIS ERROR

SAGA GIS ERROR

TOOLBOX OPTIONS SETTINGS:

Finally I think so this is solution for my another question Why can't I use processing algorithms in a standalone PyQGIS script?

qgis - Issues creating slope raster layer from a point layer

I'm trying to create a slope raster layer using QGIS 2.0 on Windows XP.

I started with a points layer with the Z elevation attribute.

I then used the interpolation plugin to create a DEM.

But, when I used the Raster Terrain Analysis plugin, it gave me a strange output. What can I be doing wrong?

postgis - How to route over (partial) edges in pgrouting?

I have made a route planning system (based on shortest path using A*) like the Google directions service using pgrouting and google maps. It includes clicking on the map to generate waypoints, dragging of waypoints or polylines to change your desired route. It also shows a directions panel with instructions how to get from A to B.

However in the reverse geocoding i use to find out where the user clicks on the map i return the closest node in the topology. This all works pretty well except for some real world cases where a direct route from A to B is longer than an indirect route through D

The blue A,D,B are the nodes and the red C-A, C-B are my mouse clicks that corresponded to those nodes. The red C-C is the place on the edge where i want my route to visit.

If a user wants to go from A to B through the long direct edge it isn't possible atm because it isn't the fastest route.

I am looking for a way to route over edges or partial edges. One can imagine that a user doesn't want to start in a node but somewhere halfway on an edge (this is another use case but i feel it will be solved as well with the answer i seek).

The cleanest way i can come up with is generating temporary nodes on an edge when a user clicks there (see click C-C on the image) by effectively splitting up an edge into 2 edges and placing a new (temporary) node on the split point. This way pgrouting can route over this new node and thus force the route to follow the edges surrounding that new node.

However, this will generate a lot of temporary nodes when dragging a polyline across the map. All those during_drag_generated nodes have to be removed except the one created at drag_end as that one becomes part of the final route.

So i am wondering are there easier ways to solve this?

Answer

I solved this problem by indeed adding a temporary node on the clicked edge and adding 2 temporary edges to this temporary node. By manually joining these temporary node and edges before calling shortest_path_astar they are being used for calculating the correct path but they dont clutter your database with extra records that are only interesting for 1 particular user (at 1 particular moment in time).

Here the SQL query i used (in Python):

"SELECT 
    row_number() over (range unbounded preceding) as rownumber, 
    astar.vertex_id, 
    astar.edge_id,
    astar.cost
FROM 
    shortest_path_astar('SELECT 
                            gid as id, 
                            source::integer as source, 

                            target::integer as target,
                            length::double precision as cost, 
                            ST_X(ST_Startpoint(the_geom)) as x1, 
                            ST_Y(ST_Startpoint(the_geom)) as y1, 
                            ST_X(ST_Endpoint(the_geom)) as x2, 
                            ST_Y(ST_Endpoint(the_geom)) as y2 
                        FROM ways 
                        %s', 
                        %s, 
                        %s, 

                        false, 
                        false) astar" % (extra_edges, route_origin, route_destination)

where route_origin or route_destination can be the new temporary (negative) id of the virtual node and extra_edges looks like:

" UNION SELECT %d, %s, %d, ww.length * %f, ST_X(ST_Startpoint(ww.the_geom)), ST_Y(ST_Startpoint(ww.the_geom)), ST_X(ST_Line_Interpolate_Point(ww.the_geom, %f)), ST_Y(ST_Line_Interpolate_Point(ww.the_geom, %f)) FROM ways ww WHERE ww.source = %s AND ww.target = %s" % (edgeid, edge_fromnode, nodeid, perc, perc, perc, edge_fromnode, edge_tonode)
+ " UNION SELECT %d, %d, %s, ww.length * (1-%f), ST_X(ST_Line_Interpolate_Point(ww.the_geom, %f)), ST_Y(ST_Line_Interpolate_Point(ww.the_geom, %f)), ST_X(ST_Endpoint(ww.the_geom)), ST_Y(ST_Endpoint(ww.the_geom)) FROM ways ww WHERE ww.source = %s AND ww.target = %s" % (edgeid-1, nodeid, edge_tonode, perc, perc, perc, edge_fromnode, edge_tonode);

Where edgeid is a unique (negative) id for this temp edge, nodeid a unique (negative) id for this virtual node, edge_fromnode the from node of the original edge, edge_tonode the to node of the original edge, perc is the percentage on the original edge (edge_fromnode, edge_tonode) of the new virtual node.

A week ago (9-1-2012) the pgrouting newsletter also touched this subject. See: Pgrouting-users Digest, Vol 40, Issue 2 -> http://lists.osgeo.org/pipermail/pgrouting-users/2012-January/000927.html

arcpy - ArcGIS Geoprocessing Script runs fine in Desktop but crashes as Geoprocessing Service?

I have a small and simple geoprocessing script that I am trying to expose as a geoprocessing service in ArcGIS Server 10 SP4. The script is below. Pardon the hard-coding...

With the hard-coded variables, this script runs perfectly fine in desktop. However, once published as gp service, it crashes in a bad way, like I can't even trap the error. The exact message in the GIS Server log - "Container process 6436 has crashed on machine westchamp24-pc."

A couple of potential causes I can eliminate right off the bat:

1. Permissions: The entire D: drive is wide open to the 'Everyone' user


2. 
   

   Licensing: ArcGIS Server is licensed at Enterprise level

   
   

   # Check out the ArcGIS Spatial Analyst extension license
   arcpy.CheckOutExtension("Spatial")
   
   
   # Execute CostDistance
   outCostDistance = CostDistance("origin",SOURCE,"#","D:/RasterStuff/backlinkRaster.tif")
   
   # Execute CostPath
   outCostPath = CostPath("destination", outCostDistance,"D:/RasterStuff/backlinkRaster.tif")
   
   # Convert Result to Polyline
   arcpy.RasterToPolyline_conversion(outCostPath, "D:/RasterStuff/costPath.shp")
   featSet = arcpy.FeatureSet()
   
   featSet.load("in_memory/costPath")
   
   arcpy.SetParameter(0,featSet)
   

   
   

   except: # do stuff with error print finally: print

   
   

Answer

I found that when I reverting back to arcgisscripting (instead of arcpy), the gp service works fine...

My code in this script is pretty vanilla - scraped almost entirely from ESRI site, so I have no idea what could be causing the problem.

python - How to extract specific information from GRIB files?

How to extract fetched GRIB data with Python and gribapi package (link goes to official ECMWF Wiki)? I tried to follow few examples from their "documentation" and do it myself, but I just cannot figure out how to retrieve only specific parameter (e.g. surface temperature or wind) for given latitude and longitude. I don't need to plot the data or visualizations, I just need to pass few parameters when reading the file and get numbers back.

So far, I tried to download original GRIB2 file from www.ftp.ncep.noaa.gov and run few Python API example scripts from mentioned ECMWF Wiki page. Since Python API documentation doesn't exist, I have some hard time to understand which methods should I use in order to:

• read the file,


• select only specific "message" (e.g. surface temperature only) if there's more "messages" in a file,



• print out the value for a specific location (using grib_find_nearest method).

Example code which I have so far, mostly copied from ECMWF examples and this answer from Getting time dimension from GRIB file?:

import traceback
import sys

# hopefully the only external package I need
from gribapi import *

# file with all possible 'messages'

INPUT = '2015121406_gfs.t06z.pgrb2.0p25.f000'
VERBOSE=1 # verbose error reporting

def example():
    f = open(INPUT)

    # location we are interested in
    lat, lon = 64.1353, -21.8952

    while 1:

        # STEP 1:
        # open downloaded GRIB2 file
        gid = grib_new_from_file(f)
        if gid is None: break

        # define the iterator (which is throughout the program the same?)
        iterid = grib_iterator_new(gid,0)

        # get the result for the nearest location
        nearest = grib_find_nearest(gid, lat, lon)[0]


        while 1:
            # STEP 2:
            # ???
            #
            # the loop goes through the whole file
            # instead just selecting messages we need beforehand...
            # fictional function which selects
            # 'TMP' (temperature) on 'surface' messages only
            # without need to iterate all of them:

            # 
            # result = grib_select_specific_message(gid, 'TMP', 'surface')

            result = grib_iterator_next(iterid)
            if not result: break

        # STEP 3:
        # profit!
        # result variable returns numbers which I process in any way I need, yay!


        # more undocumented stuff,
        # put here just because examples do it the same way
        grib_iterator_delete(iterid)
        grib_release(gid)

    f.close()

# main program function
def main():
    try:

        example()
    except GribInternalError,err:
        if VERBOSE:
            traceback.print_exc(file=sys.stderr)
        else:
            print >>sys.stderr,err.msg

        return 1

# run the program

if __name__ == "__main__":
    sys.exit(main())

The problem with this code is that it iterates trough the whole file, all the messages and all possible coordinates (whole planet, 0.25 degrees precision). If file is 300+ MB it takes a while to read it. I feel like using "iterator" here (whatever that is) is a bit wrong since I need to "cherry-pick" only specific information (few weather parameters for a particular location only) and grib_find_nearest() function doesn't need an iterator at all to select only part of the data.
This could be half-solved with filtering the data I need before downloading the file (results in much smaller file) but I'm still not sure how to do it (part of Downloading GRIB GFS files with specific filters?), but still I'd like to figure this out since it might happen that I will occasionally have to download "full" files.

Answer

I eventually ended up with ditching gribapi and switching to Met Office's Iris Python package which solves my problem in a very elegant way. Although it was a bit of a pain to install it (dependencies could be really tricky) at least the documentation is good and it is really easy to use it.

Drawing distances from one point to another using QGIS?

I want to be able to generate a distance line between two points on the map, which also has the distance as a label (similar to the second image seen on Google maps or in the first image from the red point to Aberdeen).I have used the measure line tool and it has provided myself with a distance, however I am yet to figure out how to produce a line with the classified distance labeled on it. does anyone know how to solve this issue. I am using QGIS versions 3.2

I want to record the distance from the red point to the location of Aberdeen so that when I go into {print composer} I can make a map which has a line running between the two distance with its distance e.g 100km recorded I don't have the distance recorded in my {attribute table} would I need to download that data then once collected, I would be able to draw the line?

Answer

This is pretty simple so I'll give you the steps but if you are brand new to GIS then you will need to look up more for each individual step but Google is your friend.

Create a new shapefile layer and set it's geometry to 'Line'.

Select the new line layer in the table of contents and toggle editing.

Enable snapping.

Click Add Line Feature and digitise a line from point A to B.

Open the layer's attribute table.

Open the Field Calculator.

Create a new field called 'Distance', make it a decimal number and give it an appropriate length and precision.

In the expression window type '$length'. This will give the field the length of the line in metres (depending on your coordinate reference system). To make this kilometres change the expression to '$length / 1000'.

Under the layer's properties enable editing using the 'Distance' field.

This will now display in the layout view.

openlayers 2 - Zooming the map after filter

I have wms layer and after i do cql filter on this. And now i want to see result. I want map zoomed to this result. My code.

ex=layer.getExtent();
map.zoomToExtent(ex);

But getExtent() returns wrong bounds. In Geoserver i have layer settings

Native SRS=70066
Declared SRS=900913
Reproject native to declared

In OL when create map i write projection: EPSG 900913;

What im doing wrong? I ll try do transformation

ex=layer.getExtent().transformation(map.projection, map.displayProjection);

projection=70066 displayProjection=900913

But it nothig to change.

Maybe im gonna do it whis Capabilities?

pyqgis - Switch QGIS proxy settings programatically

I am using QGIS at two offices with different proxy-servers.

At the moment I have to change the proxy settings manually every time i switch the office.

I was wondering if its possible to change the proxy-settings programatically with PyQGIS? then I could write a plugin to switch between the proxy-settings.

EDIT1:

In the meantime i found a way to change the proxy settings of QGIS but still it's not working.

With this code I can change the settings:

from PyQt4.QtCore import QUrl, QSettings
from PyQt4.QtNetwork import QNetworkRequest, QNetworkProxy
from qgis.core import QgsNetworkAccessManager

my_settings={"Proxy enabled": u'proxy/proxyEnabled', "Proxy Host ": u'proxy/proxyHost', "Proxy Port": u'proxy/proxyPort'}
fiddler={"Proxy enabled": True, "Proxy Host ": "localhost", "Proxy Port": 8888}

freiburg={"Proxy enabled": True, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}
aus={"Proxy enabled": False, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}

current_choice=aus

s = QSettings() #getting proxy from qgis options settings

for key, val in my_settings.iteritems():
    #print str(key)+":"+str(val)
    settings_key=key

    #print str(settings_key)
    # Get user defined current setting
    for key2, val2 in current_choice.iteritems():
        if key2==settings_key:
            #print key
            #print val
            settings_val=val2
    current_setting = s.value(str(val).decode('unicode-escape'))
    #print str(val).decode('unicode-escape')
    #print str(key)+": "+str(current_setting)

    s.setValue(unicode(str(val)), settings_val)
s.sync()

# procedure to set proxy if needed

proxyEnabled = s.value("proxy/proxyEnabled", "")
proxyType = s.value("proxy/proxyType", "" )
proxyHost = s.value("proxy/proxyHost", "" )
proxyPort = s.value("proxy/proxyPort", "" )
proxyUser = s.value("proxy/proxyUser", "" )

proxyPassword = s.value("proxy/proxyPassword", "" )
if proxyEnabled == "true": # test if there are proxy settings
   proxy = QNetworkProxy()
   if proxyType == "DefaultProxy":
       proxy.setType(QNetworkProxy.DefaultProxy)
   elif proxyType == "Socks5Proxy":
       proxy.setType(QNetworkProxy.Socks5Proxy)
   elif proxyType == "HttpProxy":
       proxy.setType(QNetworkProxy.HttpProxy)
   elif proxyType == "HttpCachingProxy":

       proxy.setType(QNetworkProxy.HttpCachingProxy)
   elif proxyType == "FtpCachingProxy":
       proxy.setType(QNetworkProxy.FtpCachingProxy)
   proxy.setHostName(proxyHost)
   proxy.setPort(int(proxyPort))
   proxy.setUser(proxyUser)
   proxy.setPassword(proxyPassword)
   QNetworkProxy.setApplicationProxy(proxy)

This works so far as I can see the changed settings in the QGIS UI (settings->options).

The settings are also written to the windows registry but the changes won't have any effect until i click the OK button in the QGIS settings dialog.

You can test this by setting the proxy programmatically to some proxy-settings that should prevent QGIS from accessing the internet (e.g. localhost:98765) and try to load and pan through a wms-layer. Any idea what's missing?

Edit2: I just piped the output from qgis to a file and had a look at what is going on when I change the proxy-settings using the GUI:

src/core/qgsnetworkaccessmanager.cpp: 364: (setupDefaultProxyAndCache) [9134ms] setting proxy 3 192.168.95.165:8080 /
src/core/qgsnetworkaccessmanager.cpp: 167: (setFallbackProxyAndExcludes) [0ms] proxy settings: (type:HttpProxy host: 192.168.95.165:8080, user:, password:not set

So there are two functions called (setupDefaultProxyAndCache and setFallbackProxyAndExcludes). Perhaps something like that has to be done when using pyQGIS to change the settings?

Answer

Problem solved. QgsNetworkAccessManager had to be set to the new values:

from qgis.core import *

from PyQt4.QtCore import *
from PyQt4.QtNetwork import QNetworkRequest, QNetworkProxy

my_settings={"Proxy enabled": u'proxy/proxyEnabled', "Proxy Host ": u'proxy/proxyHost', "Proxy Port": u'proxy/proxyPort'}
fiddler={"Proxy enabled": True, "Proxy Host ": "localhost", "Proxy Port": 8888}
freiburg={"Proxy enabled": True, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}
aus={"Proxy enabled": False, "Proxy Host ": "192.168.95.165", "Proxy Port": 8080}

current_choice=aus


s = QSettings() #getting proxy from qgis options settings

for key, val in my_settings.iteritems():
    settings_key=key
    for key2, val2 in current_choice.iteritems():
        if key2==settings_key:
            settings_val=val2
    current_setting = s.value(str(val).decode('unicode-escape'))
    s.setValue(unicode(str(val)), settings_val)
s.sync()


proxyEnabled = s.value("proxy/proxyEnabled", "")
proxyType = s.value("proxy/proxyType", "" )
proxyHost = s.value("proxy/proxyHost", "" )
proxyPort = s.value("proxy/proxyPort", "" )
proxyUser = s.value("proxy/proxyUser", "" )
proxyPassword = s.value("proxy/proxyPassword", "" )
proxy = QNetworkProxy()
#setting HttpPtoxy
proxy.setType(QNetworkProxy.HttpProxy)


proxy.setHostName(proxyHost)
proxy.setPort(int(proxyPort))
proxy.setUser(proxyUser)
proxy.setPassword(proxyPassword)
QNetworkProxy.setApplicationProxy(proxy)
network_manager=QgsNetworkAccessManager.instance()
stringlist= ""
network_manager.setupDefaultProxyAndCache ()
network_manager.setFallbackProxyAndExcludes(proxy, stringlist)

Using Spatial References with ArcPy after *.prj files no longer installed by ArcGIS 10.1 for Desktop?

I'm writing a python script that uses the CreateFeatureclass_management tool, without a template. One of the required parameters is a spatial reference. I used to be able to point to the system *.prj file in 9.3, but that's no longer an option, and I don't want to point to an existing shapefile's *.prj because that will cause issues when I share the script.

Now that Esri doesn't install all the *.prj files with ArcGIS, how can I set the projection in a script without referencing an existing shapefile?

Answer

As stated in the help simply choose your spatial reference by name or factory code.

# Using spatial reference name

sr = arcpy.SpatialReference("Hawaii Albers Equal Area Conic")

# Using factory code
sr = arcpy.SpatialReference(32145)

You can find the list of spatial reference factory code here:

clustering - Creating polyline-based "heatmap" from GPS tracks?

This winter I am planning to track my downhill skiing/snowboarding using a GPS. Most of my riding will occur at the same resort. I would like to be able to create a sort of "heatmap" that shows the amount of runs that I have made in a given area. As I add more and more GPS traces to my database, my goal would be to see a sort of linear heatmap of the most traveled areas. Given the nature of downhill skiing, you would expect the uphill chairlift lines will be the "hottest" because it will be the only places visited over-and-over again.

Given that 1) my track will not be the same every time and 2) the area covered by following one "run" may be a few hundred feet wide, what might be the best way to analyze this "linear" data to create a sort of heatmap? My thoughts were to buffer the lines, then intersect the polys to get a sort of Venn Diagram thing going. My preference is to use open-source technologies. I've got QGIS and PostGIS loaded and available.

UPDATE: In regards to @blah238's response, I was thinking of something that might be able to "collect" the number of passes ("runs") through an area, and then symbolize by the count. Conceptually, this would be similar to ArcGIS "Collect Events" (but for lines, not points) or Collapse Dual Lines To Centerline (but for multiple line in roughly the same area).

A more visual example of a similar concept might be a traffic-flow map, where highly-congested areas would equate to "highly-traveled" ski runs/areas:

I've read the following questions that may give some ideas, but they don't really address what I am trying to accomplish:

Clustering Trajectories (GPS data of (x,y) points) and Mining the data

Managing error with GPS routes (theoretical framework?)

WMS not loading in QGIS 2.6.1 Brighton

I'm using QGIS 2.6.1 Brighton, trying to load the following WMS layer: http://services.nconemap.com/arcgis/services/Imagery/Orthoimagery_Latest/ImageServer/WMSServer

I can successfully load the WMS in ArcGIS, and have spoken with the service providers, who have recently upgraded the service to work with the latest standard for WMS. I was advised to make sure the DPI mode is turned off, but it seems to have no effect.

When I try to add a new WMS layer, the layer appears in the browser window, when I expand the layer, it processes for a few seconds and then stops. Nothing displays in the expanded view and the layer will not add.

What do I need to do to get this layer to load? Does anyone else have success loading this layer in QGIS 2.6.1?

Answer

Works for me in 2.8. Try the following options in the connection:

qgis - How to see what lines have changed in an updated Shapefile?

I'm using my county's GIS data to update the roads in OpenStreetMap. I have a copy they published last year and one that was just published, and I'd like to find all the LineStrings that have either had their attributes or their geometry changed. This will assist me in making sure all the new and modified streets are updated in OpenStreetMap.

I'd like to do this using FOSS software, such as QGIS or Python/OGR. The street segments should have a unique identifier, so my only thought is to write a Python script which opens both Shapefiles, finds any segments in the new that aren't present in the old (added segments), and the reverse (removed segments), and then loop over matched pairs comparing their constituent coordinates to see if anything has changed.

Is this a good approach? Is there a simpler way?

openlayers 2 - Showing Two Labels of One feature In GeoServer?

I have prepared a Web GIS Map using Geoserver-Openlayers_Postgis....For One Feature I need to show two labels such as "premise_no' & "bld_name'... The following is the xml I used for styling(in geoserver).... But it is showing either/or labels...

xml for layer styling(in geoserver)....

    
        
            
        
        
            v_building_ar
            
            
                group 0

                Feature
                generic:geometry
                simple
                
                    1
                    1.0
                    2000.0
                    
                        
                            #FF8000

                            round
                            0.05
                        
                    
                    
                        
                            prem_no
                        
                        
                            Arial

                            8.0
                            normal
                            normal
                        
                        
                            
                                
                                    0.0
                                    0.0
                                

                                
                                    0.0
                                    0.0
                                
                            
                        
                        
                            #000000
                        
                    

                  
                        
                            bldname
                        
                        
                            Arial
                            8.0
                            normal
                            normal
                        

                        
                            
                                
                                    0.0
                                    0.0
                                
                                
                                    0.0
                                    0.5
                                

                            
                        
                        
                            #000000
                        
                                      
                
            
        
    

Answer

The label placement algorithm will almost certainly decide not to show both labels as they will overlap. There are various vendor options that you can use to make it less picky about overlaps and clutter but it will end up looking ugly.

The better solution is to combine the two properties into one label (this will give you much more control over the way they line up etc. too). See the user guide for full details but try something like:

   prem_no, bldname

postgis - QGIS installation fails under Ubuntu 16.04

I decided yesterday to upgrate QGIS from 2.16 to 2.18. This turned out to be a stupid idea, because the following happened:

After removing, purging and autoremoving qgis, python-gis and grass I tried to install the new version using the http://qgis.org/debian repositories.

The result was:

The following packages have unmet dependencies.
python-qgis : Depends: python-qt4-sql but it is not installable
               Depends: python-qgis-common (= 1:2.18.1+24xenial) but it is not going to be installed

               Depends: python-psycopg2 but it is not installable
               Depends: python-qscintilla2 but it is not installable
               Depends: python-jinja2 but it is not installable
               Depends: python-markupsafe but it is not installable
               Depends: python-dateutil but it is not installable
               Depends: python-requests but it is not installable
               Depends: python-tz but it is not installable
               Depends: python-yaml but it is not installable
               Depends: python-future but it is not installable
               Depends: python-pyspatialite but it is not installable

               Depends: libqgispython2.18.1 but it is not going to be installed
               Depends: libqgis-analysis2.18.1 but it is not going to be installed
               Depends: libqgis-core2.18.1 but it is not going to be installed
               Depends: libqgis-gui2.18.1 but it is not going to be installed
               Depends: libqgis-networkanalysis2.18.1 but it is not going to be installed
               Depends: libqgis-server2.18.1 but it is not going to be installed
               Depends: libqscintilla2-12v5 (>= 2.8.4) but it is not installable
               Recommends: liblwgeom-dev but it is not installable
 qgis : Depends: libgdal.so.1-1.11.3 but it is not installable
        Depends: libgdal1i (>= 1.8.0) but it is not installable

        Depends: libgeos-c1v5 (>= 3.4.2) but it is not installable
        Depends: libqgis-analysis2.18.1 but it is not going to be installed
        Depends: libqgis-app2.18.1 but it is not going to be installed
        Depends: libqgis-core2.18.1 but it is not going to be installed
        Depends: libqgis-gui2.18.1 but it is not going to be installed
        Depends: libqgis-networkanalysis2.18.1 but it is not going to be installed
        Depends: qgis-providers (= 1:2.18.1+24xenial) but it is not going to be installed
        Depends: qgis-common (= 1:2.18.1+24xenial) but it is not going to be installed
        Recommends: qgis-provider-grass but it is not going to be installed
 qgis-plugin-grass : Depends: qgis-provider-grass (= 1:2.18.1+24xenial) but it is not going to be installed

                     Depends: libgdal1i (>= 1.8.0) but it is not installable
                     Depends: libqgis-app2.18.1 but it is not going to be installed
                     Depends: libqgis-core2.18.1 but it is not going to be installed
                     Depends: libqgis-gui2.18.1 but it is not going to be installed
                     Depends: libqgisgrass7-2.18.1 but it is not going to be installed
                     Depends: grass-core but it is not installable
E: Unable to correct problems, you have held broken packages.

I checked online for help and found Broken packages and unmet dependency installation QGIS, PostgreSQL and postgis Ubuntu 14.04 and QGIS install on Ubuntu 14.04 fails

The problem described there is very similar to mine only that I am working on ubuntu 16.04 and the missed packages are slightly different. Unfortunately the proposed solutions doesn't work for my case:

• 
  

  Installing the named packages individually only leads to more uninstalable packages.

  
  


• 
  

  Using the ubuntugis repository results in the error message:

  
  

  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY 089EBE08314DF160

  
  

  and ignoring this again lead to no installable files.

  
  


• 
  
  

  Using sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys 073D307A618E5811

  
  

  only leads to:

  
  

  gpg: requesting key 618E5811 from hkp server keyserver.ubuntu.com gpg: keyserver timed out gpg: keyserver receive failed: keyserver error

  
  


• 
  

  apt list --installed | grep gdal dosn't find anything

  
  


• 
  

  And installing PostGIS manually also ends up in unistalable packages.

  
  

Can somebody please tell if there is anything else I can do?

I have been struggling here for more than a day.

I assume that something is messed up with the sources but I don't know how to find out what or how to repair it.

EDIT after comments

I now managed to get around the key-server errow and access the ubntugis ppa without problems.
(The solution for this was to properly install the ppa using sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable )

Trying to install qgis using the sources deb http://qgis.org/ubuntugis xenial main deb-src http://qgis.org/ubuntugis xenial main deb http://ppa.launchpad.net/ubuntugis/ubuntugis-unstable/ubuntu xenial main

again results in a number of unmet dependencies. One of them is libgdal20. Trying to install this (as recommended by ArneJ) leads to the the following unmet dependencies

libgdal20 : Depends: libarmadillo6 but it is not installable

             Depends: libcrypto++9v5 but it is not installable
             Depends: libdap17v5 but it is not installable
             Depends: libdapclient6v5 but it is not installable
             Depends: libepsilon1 (>= 0.8.1) but it is not installable
             Depends: libfreexl1 (>= 0.0.2~beta20110817) but it is not installable
             Depends: libgeos-c1v5 (>= 3.4.2) but it is not installable
             Depends: libgeotiff2 (>= 1.4.1) but it is not installable
             Depends: libhdf4-0-alt but it is not installable
             Depends: libhdf5-10 but it is not installable
             Depends: libkmlbase1 (>= 1.3.0~r864) but it is not installable

             Depends: libkmldom1 (>= 1.3.0~rc2) but it is not installable
             Depends: libkmlengine1 (>= 1.3.0~r864) but it is not installable
             Depends: libmysqlclient20 (>= 5.7.11) but it is not installable
             Depends: libnetcdf11 (>= 4.0.1) but it is not installable
             Depends: libogdi3.2 but it is not installable
             Depends: libopenjp2-7 (>= 2.0.0) but it is not installable
             Depends: libpq5 but it is not installable
             Depends: libproj9 (>= 4.8.0) but it is not installable
             Depends: libqhull7 but it is not installable
             Depends: libspatialite7 (>= 4.2.0) but it is not installable

             Depends: libxerces-c3.1 but it is not installable
             Recommends: proj-bin but it is not installable

Answer

Problem solved!!! The problem were in fact the missing software sources and How do I resolve unmet dependencies after adding a PPA? finally helped.

(I don't understand how I was able to install qgis16 earlier, but some mysteries have to remain unsolved...)

geoserver - OpenLayers 3: Cross-Origin Request Blocked: The Same Origin Policy disallows

Using OpenLayers 3, I cannot get this message to go away:

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON. This can be fixed by moving the resource to the same domain or enabling CORS.

This is the code:

// Ol3 only supports Projections "EPSG:4326" and "EPSG:3857". For every other projection you need proj4js
        proj4.defs("EPSG:2236", "+proj=tmerc +lat_0=24.33333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs");


        // Leases Layer
        var myLayer = new ol.layer.Vector({
            source: new ol.source.GeoJSON({
                projection: 'EPSG:2236',
                url: 'http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON',
                crossOrigin: null
            })
        });


// View
        var view = new ol.View({
            projection: 'EPSG:2236',
            center: [0, 0],
            zoom: 4
        });

        // Map
        var map = new ol.Map({
            target: 'map',

            renderer: 'canvas',
            layers: [myLayer],
            view: view
        });

I have tried setting the crossOrigin setting to:

crossOrigin: null
crossOrigin: 'null'
crossOrigin: 'anonymous'

I only see the zoom in/out control but the layer is not rendered.

---

I went with simon's option 3 below. I enabled CORS in GeoServer by copying the necessary jetty-servlets jar files and enabling it in the \WEB-INF\web.xml:

    cross-origin
    org.eclipse.jetty.servlets.CrossOriginFilter
    
        allowedOrigins
        *
    

    
        allowedMethods
        *
    
    
        allowedHeaders
        *
    



    cross-origin
    /*

After I did that, I tested the page again and receive the same error:

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON. This can be fixed by moving the resource to the same domain or enabling CORS.

Looks like I am still missing something. Do I have to do anything from the OpenLayers Side?

---

I ended up getting rid of Jetty and uninstalling GeoServer completely. The problem is when you install the geoserver windows installer, it installs a version of jetty that is 4 years old! (Jetty version 6.1.8) Even though I had copied the jar files for CORS, it is only supported in Jetty 7+.

I found out that you can install a WAR file. I decided to use Tomcat since that is what GeoServer is mostly tested on according to this note from the GeoServer website:

Note GeoServer has been mostly tested using Tomcat, and therefore these instructions may not work with other container applications.

These are the instructions for installing the WAR file:

http://docs.geoserver.org/stable/en/user/installation/war.html

This is a nice how-to video also:

https://www.youtube.com/watch?v=YEOA8WWWVCw

After you complete the install, you then enable CORS:

http://enable-cors.org/server_tomcat.html

spatial database - Need a vector format that is editable in QGIS and supports >10 character column names

For reasons, I need to create a new column in a vector shapefile with a name longer than 10 characters, and export to mapinfo TAB format. I know .shp's don't support >10 character column names, so I'm looking for an intermediate format that's both editable in QGIS, allows addition of columns after features have been added, and supports >10 character UPPERCASE column names. All of the vector formats seem to support some combination of these, but not all three.

Anybody have this problem in the past, and/or know of a format that supports this use case?

Answer

Any SQL-compliant spatial database should work, including PostGIS (if you happen to be running a PostGIS server) or SpatiaLite (but note quirks about behavior mentioned below).

Either one will satisfy your criteria:

1. Identifiers (including column names) can be up to 63 characters in a default PostGIS installation (and this can be increased by changing the server's NAMEDATALEN constant). I can't find any limit to column name length in SQLite docs, and blog posts seem to confirm that there aren't any.



2. Fields are editable in QGIS.


3. New columns can be added at any time.


4. Uppercase column names (UGGH!) can be used, but the behavior is different and complex. See following.

In general, SQL developers will try to avoid using case-sensitive identifiers. However, case-sensitive identifiers can be created using quoting. Without quoting, some databases (Oracle) will fold identifiers to UPPERCASE, some (PostgreSQL/PostGIS) will fold identifiers to lowercase, and some (SQLite/SpatiaLite, SQL Server) are case-preserving but not case-sensitive.

In PostGIS you can force uppercase by using quoted identifiers:

CREATE TABLE foo (
  "COLUMN_NAME_1" int
);

Note that since I never do this kind of table creation in QGIS DB Manager, I don't know whether you can force uppercase column names that way, but you can do so using other management tools or writing the SQL yourself (which you can submit using DB Manager's SQL Editor) as demonstrated above.

This will require you to use quoting in any SQL that you write, as SELECT COLUMN_NAME_1... will be folded to SELECT column_name_1... internally, and column_name_1 does not equal COLUMN_NAME_1.

For SpatiaLite, the behavior is a little strange. SQLite will preserve the case of the column as created (whether or not you quote the identifier), and will also preserve it on export, but in any SQL that you write the column name will be treated in a case-insensitive manner. Any combination of upper and lower case letters will be accepted regardless of quoting, and SQLite will treat COLUMN_NAME_1 and column_name_1 as a conflicting names and will not let you create both columns in the same table. (PostGIS, OTH, will allow names that differ only in case.)

For either PostGIS or SpatiaLite, the case will be preserved if you load the layer in QGIS and export to TAB using the Save As dialog.

More info on identifier case-sensitivity is available in this extremely informative blog post: http://www.alberton.info/dbms_identifiers_and_case_sensitivity.html

NOTE: Original answer only discussed PostGIS. Based on answer by @user30184, I have significantly changed answer to include information about SpatiaLite.

Making stops not in certain sequence in order to find shortest route with ArcGIS Network Analyst?

Is it possible to make stops not in a certain sequence in order to find the shortest route with the network analyst?

I now use the small green points as an input stops because I want the lines with a pink dot beside it in the route. Then I reorder them and don't preserve the end and startpoint. I thought this meant that the network analyst doesn't make a real sequence of the input stops and doesn't choose a real start and end point, and just finds the shortest route in between the points.

However if I look at the output of the tool (light green), it seems like the network analyst tool still uses a kind of sequence and begin and endpoint. Because in the picture it looks like the point where the blue arrow points to is a start or end point. The red line is the route I want. So the route makes a round now (at the left where the black arrow is) going past the south segment instead of going to the north and connect via the small grey link near the blue arrow (like the red route I want). The green output I now receive is not the shortest route.

I thought that one problem could be that the point near the blue arrow is assigned by the tool as start or end point and the other problem that the network analyst takes segment near the black arrow because the point in the south west (black circle) comes later in the sequence than the point in the middle (black circle). So, it gives each point 1,2,3 and so on, instead of tread them like all 0,0,0. So, I thought that this could be solved by having the stops as equal points (all zeros or something) and calculating the shortest route in between those instead of having the stops as a sequence and going from 1 to 2 and 3 and so on. But I don't know if this is right, see also last comment of @ChrisW, Using line feature as stops input in network analyst?.

geoserver - Are Dimensions/Attributes Other Than Time and Elevation Possible

I have Time Dimension Working with an Image Mosaic layer. Is it possible to add dimensions other than time/elevation? For example the files in my mosaic are named as follows where rtma or pris represents a "climate" attribute:

pristmax_20150115.tif
pristmax_20150115.tif
rtmatmax_20150116.tif
rtmatmax_20150116.tif

In my mosaic directory, I now have two regex property files (one for time dimension and one for climatesource dimension):

timeregex.properties containing regex=[0-9]{8}

climateregex.properties containing regex=^.{4}

I modified my indexer.properties property file look like so:

TimeAttribute=ingestion
ClimateAttribute=climatesource
Schema=*the_geom:Polygon,location:String,ingestion:java.util.Date,climatesource:String
PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](ingestion),StringFileNameExtractorSPI[climateregex](climatesource)

I also tried modifying postgis_rasters.properties to have a ClimateAttribute:

#-Automagically created from GeoTools-

#Wed Nov 25 13:16:00 MST 2015
Levels=0.026578191679850746,0.026578191679850746
Heterogeneous=false
TimeAttribute=ingestion
ClimateAttribute=climatesource
AbsolutePath=false
Name=postgis_rasters
TypeName=postgis_rasters
Caching=false
ExpandToRGB=false

LocationAttribute=location
SuggestedSPI=it.geosolutions.imageioimpl.plugins.tiff.TIFFImageReaderSpi
CheckAuxiliaryMetadata=false
LevelsNum=1

Geoserver created a table in Postgis as expected:

1 | 0103... |   rtmatmax_20150115.tif | 2015-01-14 17:00:00.000000 |    rtma

2 | 0103... |   pristmax_20150115.tif | 2015-01-14 17:00:00.000000 |    pris

But using &climate=rtma or &climate=pris in my url doesn't effect which raster is retrieved, whereas &time=2015-01-14 or &time=2015-01-15 selects the raster from that day. How can I get the climate attribute to function?

--- UPDATE ---

After following @rovo's answer, the new dimension appears in geoserver. Here is a screen cap.

And here is an example the postres side (note my regex was wrong so I currently don't have what I meant to have under climatesource but this is just a test):

The part I'm not yet able to get working is geoserver seems to be ignoring the CLIMATESOURCE param in my url.

In the getcapabilities I see

2016-01-01T00:00:00.000Z,2016-01-02T00:00:00.000Z...
agdd_20160101_base_fifty,agdd_20160101_base_thirtytwo.....

When using &time=2017-07-06&CLIMATESOURCE=agdd_20170206_base_thirtytwo I get a valid map back because 2017-07-07 happens to match the default for CLIMATESOURCE. When I change time to be anything else it doesn't return a map because it's still using the default "agdd_20170706_base_thirtytwo" even if I try overriding CLIMATESOURCE with a matching date.

--- UPDATE2 SOLVED ---

RoVo's answer now clarifies to use DIM_CLIMATESOURCE in the WMS request.

Answer

ClimateAttribute is an unknown keyword and will be ignored when the indexer.properties is read. The correct keyword to add additional dimensions to an Imagemosaic is AdditionalDomainAttributes.

The following should work adding a new Dimension in the Dimensions tab:

indexer.properties

TimeAttribute=time
AdditionalDomainAttributes=climatesource(climatesource)

PropertyCollectors=TimestampFileNameExtractorSPI[timeregex](time),StringFileNameExtractorSPI[climateregex](climatesource)
Schema=*the_geom:Polygon,location:String,time:java.util.Date,climatesource:String

climateregex.properties

regex=^.{4}

Note: String and Date dimensions are working properly, but due to a bug, Integer or Double Values are read as String only. That makes it not possible yet to query intervals. See this Bug Report.

Update: The bug has been fixed in 2.11.2.

WMS Request:

 &DIM_CLIMATESOURCE=...

qgis - Postgis 2.1 error after update

After a usual update progress - "apt-get upgrade" - I can not connect with qGIS to my PostGIS DB anymore.

Error in the qGIS logfile:

Your database has no working PostGIS support. 
Erroneous query: SELECT postgis_version() returned 7 [ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory]
Retrieval of postgis version failed

After asking for my version using SELECT PostGIS_full_version(); This occured:

ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory 
CONTEXT:  SQL statement "SELECT postgis_lib_version()" 

PL/pgSQL function postgis_full_version() line 19 at SQL statement

How can I fix that?

Thats my environment: xUbuntu 14.04, PostgreSQL 9.3, Postgis 2.1, QGIS 2.3,

Answer

As you can read here it could be a bug of the upgrade procedure.

Waiting for fixing I advise you to can follow the workaround:

apt-get remove postgresql-9.3-postgis-2.1 liblwgeom-2.1.3 liblwgeom-2.1.2 postgis

The actual versions can be other than in the example. Use commands like this to find the packages to remove:

dpkg -l | grep liblwgeom
dpkg -l | grep postgis

Followed by:

apt-get install postgresql-9.3-postgis-2.1

How to extend line features up to their intersection point in QGIS?

Is there a tool or plugin in QGIS to extend line features up to their intersection point?

Answer

CADDigitize now comes with a trim/extend function.

distance - Finding and using data on population density files?

I'm attempting to replicate this map for the state of Nevada in QGIS:

source: http://www.bostonfed.org/commdev/c&b/2014/spring/map-distance-to-closest-snap-outlet.htm

(A SNAP outlet is just a store that accepts food stamps.)

I already got my zip code and census tract layers, and geocoded the SNAP outlets in Nevada.

Looking at the map, it doesn't seem like they used the centroid of each census tract to determine distance. I think they might have used population density. This makes sense because the majority of the population in any given tract could be say, in a corner. So it doesn't make sense to use the centroid of a tract. Assuming this is correct, what I want to do is import a population density point layer. Then, by census tract, find the centroid of the population density, and determine distance to the nearest SNAP outlet from there. I was able to find a population density grid, but I don't know how to use it. Upon importing it, it blackened my map.

I came across some population data on census.gov, but it appears to be population count by area (polygon), not point data.

qgis - Raster generalization - buffers in rasters, expand pixels?

Is it possible to create buffers for pixels of raster files? Actually I need to expand the pixels classified as 1 (in white in the figure) with a spatial range of 1 or 2 pixels, in order to perform some generalization. I'm using QGIS, is this possible with gdal or something?

Answer

QGIS provides an interface to GRASS GIS, which started life as a raster GIS and therefore should provide some efficient tools to tackle this problem. Referring to its manual pages of raster commands we can find the following solutions:

r.buffer -- direct buffering of white cells.

r.cost -- can compute distances to white cells. Follow this with a comparison to select short-distance cells.

r.grow -- a local morphological operation designed specifically to expand white cells into their immediate neighbors.

r.mfilter -- a general focal filter. Various focal statistics, such as max, mean, sum, median, and standard deviation can detect the presence of white cells within local neighborhoods. Follow this with a comparison to select such cells.

r.neighbors -- an even more general focal filter, which can be used similarly to r.mfilter.

r.resample -- resampling onto a coarser grid is one way to expand the white cells. The result will be somewhat "blocky".

r.spread -- letting white cells "spread" into their neighborhoods will achieve the desired buffering.

We should expect r.buffer, r.grow, and perhaps r.mfilter to use the most efficient code. (I haven't tested these to find out.)

wmts - Easy way to see the requests QGIS makes to a remote server?

I'm attempting to debug an issue with a remote ArcGis Server's WMTS as seen from GeoServer. I can see the requests GeoTools is making (because they throw an error and get logged). However the service works in QGIS (3.0 if that makes a difference) and I would like to know what the difference in the requests is.

I know I can set up wireshark or some such but I hoped there was a way of getting QGIS to log the URLs being requested directly via a flag or setting?.

Answer

It seems there is no easy way to do this in QGis, so I solved it by adding a simple proxy, Python Logging Proxy which "does what it says on the can".

I modified the file LoggingProxyHTTPHandler.py to comment out line 105

print response.content

since I didn't need to see the "contents" of the images being sent back. Then to run the proxy all you do is

python ./proxy.py

And in QGis go to Settigns->Options->Network and turn on the use proxy for web access and fill in localhost and 8000 for the Host and Port.

Then every time QGis requests a WMTS (or other web image) you will see the request and the response in the terminal.

Is there any other option to display DMS grid and graticule in QGIS Print Composer?

Is there any other option to display DMS grid and graticule in QGIS Print Composer, since there is no plugins available yet?

arcpy - Create a line perpendicular to an existing line in ArcGIS

For an existing set of line features, I would like to generate a perpendicular line, of say 1km, at the midpoint of each existing line feature in ArcGIS.

I imagine that The methodology would be to:

1. create a tangent line at the midpoint of each existing line feature


2. calculate the angle of the tangent line


3. based on the angle of the tangent line, calculate the location of the endpoint of the new line


4. create the new line with the calculated endpoint of the new line and the existing midpoint of the old line

Let's assume that we are in a UTM projection and can use coordinates to calculate distance and direction.

How would I do this with ArcGIS tools or with ArcPy?

Answer

• Create midpoints, using one of multiple possible techniques.


• Buffer them by small number, e.g. 0.5 m


• Clip originals by buffer, output - SHAPEFILE

Use this field calculator expression on field Shape:

def RotateExtend(plyP,sLength):

 l=plyP.length
 ptX=plyP.positionAlongLine (l/2).firstPoint
 ptX0=plyP.firstPoint
 ptX1=plyP.lastPoint
 dX=float(ptX1.X)-float(ptX0.X)
 dY=float(ptX1.Y)-float(ptX0.Y)
 lenV=math.sqrt(dX*dX+dY*dY)
 sX=-dY*sLength/lenV;sY=dX*sLength/lenV
 leftP=arcpy.Point(ptX.X+sX,ptX.Y+sY)
 rightP=arcpy.Point(ptX.X-sX, ptX.Y-sY)

 array = arcpy.Array([leftP,rightP])
 section=arcpy.Polyline(array)
 return section

RotateExtend( !Shape!, 1000)

UPDATE:

Ooops! it creates perpendicular, not tangent. Can be easily modified

Cleaning/editing GPX (tracks) files recorded from GPS device in QGIS?

During a long trip, I have recorded tracks of my journey with the GPS of my smartphone. I have ended up with about 30 GPX datafiles.

Now that I am back, I am trying to draw a map of this trip. I imported each file in QGIS as vector layer, and visualized the tracks with OpenLayers Plugin.

Now I need to edit some of these files:

• basically trimming the end and beginning of many files


• merging some tracks


• adding a title or a description for each track

What is the way to do it with QGIS ?

I would like to keep as much data as possible from the source data (timestamp, elevation)

Answer

The best solution I have found so far is : QLandkarte GT

It does all the actions I have listed and much more. The only trick is to get map overlay, but it is easier with the latest version.

It is worth to mention that it's multiplateform software.

Setting composer layout units to inches in QGIS

I open a new composer and select ANSI A for my paper. The notation after the name is "Letter; 8.5inx11in", but the width and height are displayed in mm, and the units selection is grayed out, so I can't change it. The rulers at the top and left of the layout are in mm, and the parameters for "Snap to grid" are also in mm. If I then select "Custom", I can select inches for units, and the Width and Height display as 8.5 and 11. However the rulers in the layout are still in mm, as are the units for the grid.

Can the entire composer project be converted to inches?

qgis 2.01 64 bit
windows 7, 64 bit

qgis - Seeking method for automatic aggregation of polys with sparse data?

I am part of a project on presenting health related data on maps. In order to maximize the analytical usefulness of the data, we wish to divide them into as small areas as possible. However, to protect individuals we also have to make sure that each area contains enough cases to secure anonymity. We are therefore looking for methods to automatically aggregate an area with too sparse data to a neighboring area – or perhaps a bit more advanced – to a non-neighboring area in the map, which according to sociodemographic parameters resembles the area in question.

I seek any input on methods for automatically:

1. 
   

   Aggregate an area with too few cases to a neighboring area – for example the neighboring area with the lowest number of cases.

   
   


2. 
   

   Aggregate an area A to the area B – neighboring or non-neighboring – which on a number of parameters has the closest resemblance to area A.

   
   

The work will be done in QGIS and/or ArcGIS Desktop, so any readily implemented methods in any of these systems are welcome, but general theoretical descriptions of statistical methods – probably most relevant for the second case – are also of interest.

Answer

Further research came across this question, which terms what you are trying to do not only as aggregation, but more specifically Zone Design. There may already be tools out there designed to aggregate data for solving this specific type of problem. Based on an answer in the aforelinked question, AZTool may fit your needs. If not, this will require a model or script to completely solve 'automatically'.

I agree with GISKid's comment that merging to non-neighboring areas would be a mistake given you are looking at a geographic relationship component. No matter how many other parameters match, if they're in completely different places you're changing distribution rather than aggregating. Much beyond that I can't address the statistical rationale for making certain choices; only the way to carry them out.

The specific attributes and structure of your data is going to play a significant role here, so I'm just going to outline a general process/workflow. I'm assuming you have a single polygon layer with attributes.

1. Create a new attribute field called ToMerge with a default value of 0.



2. Select all records where cases < cut-off, and set ToMerge to 1.


3. Iterate through each of those records to find a poly to merge to. This is the tricky part, may actually require multiple steps or a significant decision tree, and you'll have to make the rules to. If two are already adjacent, just merge them? Three or more adjacent? Merge to the neighbor with the lowest count? The highest? Match specific parameters? Check for more than one condition? Iterate through each of the neighbors and find the best match, or take the first one? Account for polygon size (significant alteration of density - say a very large poly with too few cases next to a very small poly with a lot of cases)? Preserve case density? However you decide, eventually you'll pick one neighbor and change its ToMerge value to 1 before moving on to the next feature that needs a neighbor selected.


4. After going through all of step 2's records and assigning a neighbor to be merged to, run a Dissolve tool on the layer using ToMerge as the dissolve field, do not allow multipart features, and set up the required statistics fields so that your attributes are summed or whatever method is required to combine the two poly's attribute values. Note I have ArcGIS's Dissolve tool in mind as I write this.


5. With the resulting dataset, ensure that there is still a ToMerge field and reset all values to 0, then repeat from step 2. The repeat may not be necessary if you take care of the check in step 3, but you do need to make sure that your aggregated/dissolved polygons meet the cut-off even after adding two together.


6. Once all loops are complete, you'll have your output dataset that has been aggregated as required.

Another approach rather than Dissolve (again ArcGIS) would be the Aggregate Polygons tool, followed by some joins and field calculations to aggregate the attribute values. You'd still have to have pick which ones to aggregate though.