texts: Can't get my QGIS external map to display on Ubuntu trusty: QObject::connect: Cannot connect (null) error

texts:
I've been following along the examples from the book "Building Mapping Applications with QGIS and for the life of me I can't figure out why my shapefiles aren't loading and displaying in the app window under a Ubuntu system.  The same code runs fine with my Macbook Pro setup.

I use a run.sh file to set up the linux environment:

export PYTHONPATH="/usr/share/qgis/python"
export LD_LIBRARY_path="/usr/share/qgis/lib"
export QGIS_PREFIX="/usr/share/qgis"
python externalApp.py

I have tested that I can import all qgis, qgis.core, qgis.gui etc files without errors in Python.

The Linux script I'm using is pasted below. I have played around with the shapefile name and location.

import sys
import os
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import Qt

#############################################################################

class MapViewer(QMainWindow):
    def __init__(self, shapefile):
        QMainWindow.__init__(self)
        self.setWindowTitle("Map Viewer")

        canvas = QgsMapCanvas()
        canvas.useImageToRender(False)
        canvas.setCanvasColor(Qt.white)
        canvas.show()

        layer = QgsVectorLayer(shapefile, "layer1", "ogr")
        #if not layer.isValid():
        #    raise IOError("Invalid shapefile")

        QgsMapLayerRegistry.instance().addMapLayer(layer)
        canvas.setExtent(layer.extent())
        canvas.setLayerSet([QgsMapCanvasLayer(layer)])

        layout = QVBoxLayout()
        layout.addWidget(canvas)

        contents = QWidget()
        contents.setLayout(layout)
        self.setCentralWidget(contents)

#############################################################################

def main():
    """  Our main program.
    """
    QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'], True)
    QgsApplication.initQgis()

    app = QApplication(sys.argv)

    viewer = MapViewer("/home/cary/workspace/temp/alaska.shp")
    viewer.show()

    app.exec_()

    QgsApplication.exitQgis()

#############################################################################

if __name__ == "__main__":
    main()

Computers used:
1) Toshiba Satellite windows machine, reconfigured with ubuntu 14.04
2) Macbook Pro computer

The Toshiba Satellite machine was a new machine and newly configured with Ubuntu.  I basically just installed Ubuntu then updated the system and then installed QGIS with this script:

sudo sh -c 'echo "deb http://qgis.org/debian trusty main" >> /etc/apt/sources.list'

sudo sh -c 'echo "deb-src http://qgis.org/debian trusty main" >> /etc/apt/sources.list'

gpg --keyserver keyserver.ubuntu.com --recv 47765B75

gpg --export --armor 47765B75 | sudo apt-key add -

sudo apt-get update && sudo apt-get install qgis python-qgis

Setup for both computers:
Both computers are using Python 2.7.6

Standard setups for both of the machines.  I pretty much just ran the script above for the Linux box.

Verified that the "/usr/share/qgis" folders were setup in Linux properly
Verified that the structures on the Mac all matched up with the description of folders and contents from the book on page 21 for setting up PYTHONPATH, QGIS_PREFIX etc...

The Linux file structures were different:
in the book:  PYTHONPATH="/usr/share/qgis/build/output/python/" while I noticed on my system that the "/build/output/" directories were missing.

The same code worked on the Mac setup but not the Linux setup.
I also verified that the path to the shapefiles were correct for the linux system.

Fortunately the Mac setup works and I can continue with the rest of the tutorial on that box but I would love to know why the Ubuntu box doesn't work.


texts: Batch export different map compositions based on the same Atlas

texts:
Assume I have an Atlas with the following elements:

'001'
'002'
'003'

I also have several compositions named:

'A'
'B'
'C'

With the print Atlas function, I can print a PDF per each composition, having all Atlas elements: that would produce something like:

A.pdf
 - 001
 - 002
 - 003

I wish, instead, to export an Atlas element with all compositions as its pages, like this:

001.pdf
 - A
 - B
 - C

Can someone point me in the right direction to programmatically do this? I've researched plugins and found Maps Printer - its great, but for automating the first process, not the second one. I've come across Python, maybe that's the way to do it?


texts: how can i set place mark with text in QGIS?

texts:
I want to mark a point in map  and place text near it. I used QgsVertexMarker to mark point .As far as now i am using tool-tip to show text. but i need clean code to show my text on map.

Following code is used to set marker at point A.

QgsVertexMarker * marker = new QgsVertexMarker(mPlugin->iface()->mapCanvas());
marker->setColor(Qt::black);
marker->setIconType(2);
marker->setCenter(A);

I need to display a text near that point . how to do by c++ code?


texts: How to overcome a map distortion when I am changing the angle in Qgis?

texts:
In this project I needed to change the map orientation in 90 degrees. I used a nice plugin "Dockable MirrorMap" that gives the opinion to change each layer separately.
The problem is that the map from the "OpenLayers Plugin" get distorted. This phenomenon also occurred when I changed the angle in the Layer Styling  bar or in the composer.
Me solution was to save the base map as a raster and then georeferencing it and change the angle. The outcome was good.
However I am wonder: is there is a way to solved this problem in a faster way?


texts: Trouble installing GDAL in R for MODIS package

texts:
I'm having trouble installing GDAL in R with the MODIS package. Every time run MODIS I get this error:

GDAL : Not available. Use 'MODIS:::checkTools('GDAL')' for more information!

I have tried using the link below, but the only thing it fixes is the MRT issue.
GDAL not available in MODIS R package?

What should I do to fix this issue.


texts: How to create buffer on one side of a stream in ArcGIS?

texts:
I need to create a buffer around waterfront park sites. But I only want "half" of the buffer --- the "half" that coincides with the side of the water body where to park is.

In addition, the buffer of one park may overlap another due to their proximity. I would like to have just one "half" buffer for each park.

I tried to buffer the parks, union them with waterbody, erase waterbody, select by location the buffer using the parks. But since the buffers overlap, they create multipart features and only one or two parts of the buffer are selected.

Below is example of what I am talking about.

For the park (green) on the left I would like a half buffer to the north (assuming north is up) and no buffer south of the waterbody.

For the parks on the right (one big and one small) I would like two separate buffers that are each single part features.

ps these are both polygon vector files and I have an advanced license

EDIT:
@ChrisW s method:


texts: How change legend icons size on leaflet map generated using QGIS2WEB plugin?

texts: I need to know how to change size of icons from the leaflet web map legend generated using qgis2web plugin. Is there any icon size information inside code which I can change to larger number? I have searched for it in the map folder, but with no success.


texts: Export TileMill project as SVG image with CartoCSS applied to it

texts:
I am able to export TileMill project as mbtiles and rendering tiles on map using leaflet which is working fine. But if I export the project as SVG or PNG image the CartoCSS which I have written is not applying for the image.
So how to export tilemill project as SVG image with CartoCSS applied to it.

This is the image in tileMill

This is the exported png image.


texts: Finding text object variable names in MapInfo/MapBasic

texts:
In MapBasic you can create a text object in a layout window by using the "Create text" statement. You can also use the same statement to update this object.

I am working on a tool that creates several plans by repeatedly exporting a layout window as an image file. Between each export I would like to update the text in the layout window to indicate the reference number of the site.

How do I go about finding the object's variable name?


texts: How to calculate the average spatial variation of GPS points in relation to the road network in QGIS

texts:
I have a vector dataset containing GPS locations/points of different vehicles and a road network dataset. Due to the inaccuracy of the GPS receiver, these gps points can deviate up to 10-20 meters from the road network (see below). To give an indication of the accuracy of the GPS signal I would like to know what the average spatial variation of these GPS points to the nearest road segment is.

Does anyone have an idea on how to calculate this in QGIS?


texts: Comparing differences in two polygon shapefiles using ArcGIS Desktop?

texts: I am a new user to Arcmap 10.  I've been trying to figure out how to compare two polygon shapefiles.  I need to figure out what the difference is from 2007 - 2013 (increase/decrease) in a couple of TMK maps.


texts:
I need to develop an interactive map (Show map and point with some related data) with OpenLayers from client side but from server side I don't know what I should use. GeoServer seems to be much more than what I need.

What I need is a service that accepts and it responds with WFS. Considering that I have access to the database that stores the transactional and georeferenced data.

Can I develop my own implementation of web service in PHP? Are there existing libraries? If not, what do you suggest?


texts: arcmap editor attribute window fill value programmatically

texts:
In my map I have a layer named PARCEL, when create new feature on this layer, the attribute window will open automatically. There are 2 attributes in this layer: PARCEL_ID and SUB_ID.

Since the PARCEL_ID is the left part of SUB_ID, I intend to display the PARCEL_ID value in SUB_ID field after user input the PARCEL_ID value, this is to minimize user key in wrong information.

Is there a way to do this without a customized editor tool?

I am using vb.net arcobjects 10.2.2.

EDIT1

The reason for use the attribute window in editor is when I tried to use update feature, it didn't work.

Actally my basic idea is after user create new feature:

(1) use max OBJECTID to identify the new created feature

(2) provided a window to ask user to input the PARCEL_ID and SUB_ID of the new created feature

(3) use feature.store to update the PARCEL_ID and SUB_ID values

But I failed to make it work. My code is like followings:

Private Sub GetMaxFID()
    Try
        Dim pMxDoc As IMxDocument, pFLayer As IFeatureLayer
        Dim pCursor As ICursor
        pMxDoc = My.ArcMap.Application.Document

        Dim pMap As IMap
        Dim pActiveView As IActiveView

        pMap = pMxDoc.FocusMap
        pActiveView = pMxDoc.FocusMap

        Dim layerNum = GetIndexNumberFromLayerName(pActiveView, My.Settings.ParcelLayer)
        pFLayer = pMxDoc.FocusMap.Layer(layerNum)
        pCursor = pFLayer.Search(Nothing, False)

        If TypeOf pCursor Is IFeatureCursor Then
            Dim pData As IDataStatistics = New DataStatisticsClass
            pData.Field = "OBJECTID"
            pData.Cursor = pCursor
            Dim pStatResults As IStatisticsResults = pData.Statistics
            newOBJECTID = pStatResults.Maximum
        End If

    Catch ex As Exception
        logger.Error(ex.ToString)
    End Try
End Sub

Protected Overrides Sub OnClick()
    Try
        Dim pMxDoc As IMxDocument
        Dim pMap As IMap
        Dim pActiveView As IActiveView

        pMxDoc = My.ArcMap.Application.Document
        pMap = pMxDoc.FocusMap
        pActiveView = pMxDoc.FocusMap

        Dim layerNum = GetIndexNumberFromLayerName(pActiveView, My.Settings.ParcelLayer)
        Dim pFLayer As IFeatureLayer

        pFLayer = pMap.Layer(layerNum)

        GetMaxFID()

        Dim featureClass As ESRI.ArcGIS.Geodatabase.IFeatureClass = pFLayer.FeatureClass

        If featureClass Is Nothing Then
            Return
        End If

        Dim feature As IFeature = pFLayer.FeatureClass.GetFeature(newOBJECTID)  ' is this line correct?

        If feature Is Nothing Then
            Return
        End If

        Dim objectIDFieldIndex As Integer = pFLayer.FeatureClass.FindField(newOBJECTID)

        feature.Value(1) = "test123"   'This line I want to set the PARCEL_ID of the new created feature.
                                       'is it correct to use .value(1) here? what should be the value of (1) in this line?
        feature.Value(2) = "test123-12" 'This line I want to set the SUB_ID of the new created feature.
                                       'is it correct to use .value(2) here? what should be the value of (1) in this line?

        feature.Store()

        'MsgBox("update finished")

    Catch ex As Exception
        logger.Error(ex.ToString)
    End Try
End Sub

I believe there is something wrong with my code, how can I make it work?


texts: Mapbox GL JS Toggle Layers default to "off"

texts:
Using Mapbox's 'Toggle Layer' example, how would I load the map with the "contours" and "museums" layers turned off?

I want the layers to be turned "off" by default and have the user turn "on" the layers.

I've tried replacing

link.className = 'active';

With

link.className = '';

Per this question, but still doesn't work.


texts: In ArcGIS Online is it possible to create a feature layer (hosted, view) from a join of a join?

texts:
In ArcGIS Online I have a feature layer with the following relationships:

parent -> child -> grandchild

I'd like to join all three tables into one (live) feature layer (hosted, view).

The first join, parent -> child, can be saved as a feature layer by checking the option "Create results as hosted feature layer view".

When I try to join this resulting feature layer to the grandchild table however, the same checkbox is greyed out.

Is what I'm trying to do even possible in ArcGIS Online?


texts: Using field from another table in popup of ArcGIS Web AppBuilder

texts:
I have an Web AppBuilder application, in popup windows it shows up to 10 fields.

For one field of it,  can I get a value from another service/SQL table by customizing application.

I tried to customize it in Visual Studio, but the solution didn't support due to HTTP/HTTPS errors.

For example actual popup window as "ID1" field from "table1" and by customizing I want to change it like:

SELECT ID2 FROM table2 WHERE ID1='@ID1'"

Get the "ID2" value from "table2" and display it in popup when opened instead of displaying "ID1".


texts: Full instruction for download OSM data in Ubuntu 12.04 LTC

texts:
Is there instruction for:

1) install software

2) download OSM data in Postgis database

3) set styles

4) generate tiles

5) initilaze MapProxy   (all in OS Linux)

6) View map (generated tiles) in ArcGis Server (or Geoserver) (in OS Windows)


texts: Mbtiles from dedicated server to mapbox account

texts:
I process Landsat Imagery on a dedicated server and generate mbtiles.
Is there a way to upload the mbtiles files straight from my server to my Mapbox account?


texts: Options for displaying 3D data using ESRI web APIs

texts:
I've been looking through the documentation of the various ESRI web APIs to see if there is any support for displaying 3D data, but haven't found anything useful. Is there any support for this - or do any other web APIs (Google, Bing, OpenLayers, etc) offer anything extra?

I'm assuming what's on offer is probably quite limited considering the need to purchase ArcGIS 3D Analyst for the desktop side of things.

I'm interesting in displaying either raster or vector data - or ideally both if that's possible. Any examples of what is possible would be extremely useful.


texts: What causes AttributeError: 'module' object has no attribute 'da'?

texts:
I have a shapefile "selected_polygons.shp" that contains 2000 polygons. I would like to select each polygon contained in the shapefile by using a search cursor in order to apply the function "CostDistance" for each polygon.

Here is my code:

import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
arcpy.env.workspace = "C:\Users\LCP"
arcpy.env.extent = "costs.tif"

with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInPatches_start:

for rowStart in rowsInPatches_start:

 ID_patch_start = rowStart.getValue("GRIDCODE")

 ## Define SQL expression for the fonction Select Layer By Attribute
 expressionForSelectInPatches_start = "GRIDCODE=%s" % (ID_patch_start)

 # Process: Select Layer By Attribute in Patches_start
 arcpy.MakeFeatureLayer_management("selected_patches.shp", "Selected_patch_start", expressionForSelectInPatches_start)

 # Process: Cost Distance
 outCostDist=CostDistance("Selected_patch_start", "costs.tif", "", "outCostLink.tif")

 # Save the output
 outCostDist.save("outCostDist.tif")

The problem is that I obtain a error message when I use:

with arcpy.da.SearchCursor("selected_polygons.shp",("*")) as rowsInPatches_start:

Here is my error message:

 AttributeError: 'module' object has no attribute 'da'


texts: Find CRS (coordinate reference system) based on known location and coordinates?

texts:
Is there any tool or any tutorial which would help me to figure out CRS (coordinate reference system) from related information?

In my particular case I know only the system is in meters, has origin (probably of zone) at E 91°30' and N 50°20' (not 100% sure about exactly 50°20') and it is presumably commonly used in Russia.

I would be glad for general solution if any.


texts: Select all polygons that are completely within a circular buffer

texts: I made several buffers around some points. Now I want to select all polygons that are completely within this buffer. My problem is that I also get the polygons that are partially outside of the buffer.


texts: How to calculate OSM WGS84 Lat and Lon to meters in C++?

texts:
Summary

I am converting OSM coordinates into X and Y values with accurate distances in meters between all points using c++. The calculations/code I am using is resulting in an incorrect ratio/scale in x and y. Where have I gone wrong in my calculation or code? I know there is an error in ratio due to roundabouts that are in reality round are showing as elliptical. I know the scale is wrong because ring roads are showing as being smaller than they should be relative to UK DSM Lidar data.

What I have tried so far is step 1. Convert (geodetic) latitude/longitude to ECEF geocentric cartesian coordinates from here https://www.movable-type.co.uk/scripts/latlong-os-gridref.html.

I have also checked the formula against this site https://gssc.esa.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion.

I have also tried applying the Helmert transform after the initial ECEF conversion to transform being the data more in line with OSGB36.

I have tried applying a different degrees to radians calculation only for my lattitude to correct the aspect ratio issue. This works to an extent in a rough way, but lacks accuracy and doesn't solve the underlying problem with my calculations or the scale.

I originally had height set as 0, but my calculation results differed from those on http://www.apsalin.com/convert-geodetic-to-cartesian.aspx, through trial an error, I adjusted the height value in my code until it produced results that matched the online converter to something like 4 or 5 decimal places.

WGS84 to Cartesian ECEF just appears (from my results) to be an orthographic projection of the WGS84 coordinates, rather than a projection that compensates for projecting onto a flat surface to give accurate distance between points.

(I have completely changed the way I have worded and presented the original question to reflect the further work I have been doing so solve my issue.)

//Constants for WGS84 to ECEF

double height = (27716.480814599*2);
double SemiMajorAxisA = 6378137.0;
double RadiusOfCurvature;
double FlatteningFactoroftheEarth = (1/298.257223563);
double EccentricitySquared = 0.00669437999014;
const double DEG_TO_RAD = 0.01745329251994330;
const double SEC_TO_RAD = 0.000004848;

//Helmerts constants for OSGB86//

double cx = -446.448;
double cy = 125.157;
double cz = -542.06;
double s = 0.0000204894;
double rx = -0.1502*SEC_TO_RAD;
double ry = -0.247*SEC_TO_RAD;
double rz = -0.8421*SEC_TO_RAD;

//WGS84 to ECEF for maximum lat and lon

MaxLatdouble = MaxLatdouble * DEG_TO_RAD;
MaxLondouble = MaxLondouble * DEG_TO_RAD;

    RadiusOfCurvature = SemiMajorAxisA / (sqrt(1 + (EccentricitySquared * (sin(MaxLatdouble) * sin(MaxLatdouble) ))));

    double MaxLatdoubleA = (RadiusOfCurvature + height) *(cos(MaxLatdouble) * cos(MaxLondouble));
    double MaxLondoubleA = (RadiusOfCurvature + height)* (cos(MaxLatdouble) * sin(MaxLondouble));

//Apply Helmert transform to max lat and lon

    MaxLatdoubleA = cx + (MaxLatdoubleA * (1 + s)) + (-rz * MaxLondoubleA) + (ry * height);
    MaxLondoubleA = cy + (rz * MaxLatdoubleA) + (MaxLondoubleA * (1 + s)) + (-rx * height);

//Now perform the same calculations used above for each OSM Lat Lon.

// .. same code as above in a loop for each coordinate

// final lat and lon for each is (maxlat - indivudual lat) and (maxlon - individual lon)


texts: Cannot Save Coordinates to Database (Mapbox API)

texts:
I am creating a crime mapping system using Mapbox API, so far I already able to manually input coordinates to the database then show them on the map. But what I want to happen is to automatically add the coordinates by just dragging or clicking the marker on the map. I tried the code below to get the coordinates wherever I drag or click the map then show it to an  field and it is already working. The only problem is whenever I click on the map and try saving it to the database the coordinate is always empty. I do not know why it is always saving as empty even the input field shows the coordinate I click on the map.
Here is the js code (Maybe there is something wrong in the code??):
<script>
var saved_markers = <?= get_saved_locations() ?>;
        var user_location = [77.216721,28.644800];
        mapboxgl.accessToken = '';
        var map = new mapboxgl.Map({
            container: 'map',
            style: 'mapbox://styles/mapbox/streets-v9',
            center: user_location,
            zoom: 5
        });

        var geocoder = new MapboxGeocoder({
            accessToken: mapboxgl.accessToken,
            //country: 'IN',
        });

        var marker ;

        map.on('load', function() {
            addMarker(user_location,'load');
            add_markers(saved_markers);
        });
        map.on('click', function (e) {
            marker.remove();
            addMarker(e.lngLat,'click');
            document.getElementById("lat").value = e.lngLat.lat;
            document.getElementById("lng").value = e.lngLat.lng;
        });

        function addMarker(ltlng,event) {
            if(event === 'click'){
                user_location = ltlng;
            }
            marker = new mapboxgl.Marker({draggable: true,color:"#d02922"})
                .setLngLat(user_location)
                .addTo(map)
                 .on('dragend', onDragEnd);
        }
        function add_markers(coordinates) {
            var geojson = (saved_markers == coordinates ? saved_markers : '');
            geojson.forEach(function (marker) {
                new mapboxgl.Marker()
                    .setLngLat(marker)
                    .addTo(map);
            });
        }

         function onDragEnd() {
            var lngLat = marker.getLngLat();
            document.getElementById("lat").value = lngLat.lat;
            document.getElementById("lng").value = lngLat.lng;
            console.log('lng: ' + lngLat.lng + '<br />lat: ' + lngLat.lat);
        }
        </script>

Thanks!


texts: How to use cartodb.js with API key for privacy of data

texts:
I would like to use cartodb.js in my web page reproduce a visualization that I created at CartoDB.  But I want to keep private my visualization and data at CartoDB (and I am going to secure my web page).

Is there a way to use the APIs in cartodb.js with my CartoDB API key somehow to reproduce in my web page the private visualization that I see on my CartoDB site?

I've been studying the documentation, examples, and source code for cartodb.js.  As far as I can tell, cartodb.js does not really provide for use of my API key, but perhaps I am missing something.


texts: Is LANDSAT 8 Resolution too low for coastline difference measurement?

texts:
I recently extracted coastlines (for my master thesis) for a Bangladesh region for the years 2013 and 2014. Regardless of the method I used, the distance between the 2 coastlines was always the same: 42 m.

This distance is way too big for a one year difference (same time of day, so tides were more or less the same). I realized that the 42 m are exactly the length of a diagonal with 30 m sides (which is exactly the resolution of these Landsat images).

In conclusion, no distances smaller than the pixel resolution can be measured. Does anyone know of any method to improve this, or free higher resolution alternatives to Landsat?


texts: OpenLayers 4 popup, prevent pointer to target data behind a div

texts:
I have a problem with the Popup of OpenLayers 4, when I open the popup info and I'm inside the div with the pointer, it is still triggering the hover popup with the data behind the popup, how can I prevent that?

I'm using this CSS code for the popup:

#popup {
  z-index: 20;
  background-color: #fff;
}

The popup is created with the overlay layer:

overlay = new ol.Overlay({
  element: container,
  autoPan: true,
  autoPanAnimation: {
    duration: 250
  }
});

But the pointer still passes through the div. the popup contains a lot of div and inside span

In the screenshot the hover popup was triggered when the pointer was on the Vehs text

You can refer to this for the code I use for the popup, the only difference is that I have a hover popup, I think the problem is on the creation of the click popup (the one with button) and the real z-index of it.

http://openlayers.org/en/latest/examples/popup.html


texts: I am very new to QGIS. I need any kind of tutorial where I can plot the overlapping areas inside the maps with different colors as attached in as an example image.


texts: Problem overlaying WFS on top of Bing maps using OpenLayers

texts:
I have a simple set of 99 polygons loaded into GeoServer from PostGIS. They were loaded from GeoJSON using ogr2ogr. (I'm using OpenGeo Suite 3.1 if it matters.) When I go to the layer preview and select OpenLayers it does display successfully. When I copy pasted the code and adjusted it slightly it also worked on the different server. However, now I'm trying to overlay it on top of a Bing base map. The base map displays just fine. The overlay does not display. (That is in Chrome, in Firebug it fails on map.addLayer() and aborts the entire script.)

Based on that walk through it would seem that the problem is somewhere in my declaration of the vector layer. BTW, I went into the geoserver workspace area and explicitly set the URI to what is listed under "featureNS." I've tried both declaring and not declaring the geometry name explicitly.

function init()
{
    var map = new OpenLayers.Map("map");

    var road = new OpenLayers.Layer.Bing({
         name: "Road",
         key: apiKey,
         type: "Road"
        });
        var aerial = new OpenLayers.Layer.Bing({
            name: "Aerial",
            key: apiKey,
            type: "Aerial"
        });//,{isBaseLayer: false});*/
        var hybrid = new OpenLayers.Layer.Bing({
            key: apiKey,
            type: "AerialWithLabels",
            name: "Bing Aerial With Labels"
        });

        map.addLayers([road, hybrid, aerial]);

    var data = new OpenLayers.Layer.Vector(
    "WFS Display",
    {
        protocol: new OpenLayers.Protocol.WFS ({
        version: "1.1.0",
        url: "http://<servername:port>/geoserver/opengeo/wfs",
        featurePrefix: "opengeo",
        featureType: "ogrgeojson1",
        featureNS: "http://<servername:port>/data/opengeo"
    }),
        strategies: [new OpenLayers.Strategy.Fixed()]
    });

    map.addLayer(data);
}

Any thoughts? Just for the heck of it I tried using OSM instead of Bing (although that's the job) and still no joy. So it isn't Bing at least. What did I do stupid?


texts: QGIS failure in PostGIS connection: "ident authentication failed for user"

texts:
I am trying to use the 'ident' authentication of for local PostGIS (2.5) server, so that each user can connect to his/her account without having to provide password.

The local connection (to 127.0.0.1) worked if I use trust in /etc/postgresql/11/main/pg_hba.conf:

host    all             all             127.0.0.1/32            trust
host    all             all             ::1/128                 trust

, but fails when trust is replaced with ident and if connection is made from QGIS (3.8.1).

ident authentication failed for user

The strange thing is that I can connect successfully (with ident authentication) without giving user names or password, both from the command line client psql and from pyscopg2 using python. But QGIS3 connection fails.

From documentation:

  ident

  Obtain the operating system user name of the client by contacting the
  ident server on the client and check if it matches the requested
  database user name. Ident authentication can only be used on TCP/IP
  connections. When specified for local connections, peer authentication
  will be used instead. See Section 19.3.5 for details.

  peer

  Obtain the client's operating system user name from the operating
  system and check if it matches the requested database user name. This
  is only available for local connections. See Section 19.3.6 for
  details.

Is this a bug in QGIS? Or am I missing something in the QGIS or PostGIS configuration?

How can I get ident working for QGIS?

Note: I googled the issue, some suggest using other authentication method like md5. But I specifically need to ident behavior for multiple users. md5 requires password, and trust is not secure. I need each user to access his/her own database without password, without allowing access to any other databases.


texts: OSM2PGSQL - Openstreetmap Tiles showing no detail at certain levels

texts:
I recently built myself a Map server using this guide:

https://www.linuxbabe.com/linux-server/openstreetmap-tile-server-ubuntu-16-04

This has worked really well to an extent.

I have however noticed peculiarities at certain zoom levels.

At this level I can see everything perfectly:

However if I zoom out further I lose all detail and am basically left with land outlines:

Is this perhaps a caching or stylesheet issue?

I'm using the carto stylesheet as suggested in the tutorial. Also if it were caching I'd expect to just see blank tiles rather than outlines. Curiously if I'm zoomed in further and the tiles haven't cached they show as blank before rendering normally.


texts: How do I include 64 bit python with ArcGIS 10.4 install?

texts:
I recently updated to 10.4 and was hoping the 64 bit python install was automatic, but apparently it isn't. I can't remember what I did in 10.3.x to get the 64 bit python. Does anyone know how to install it so that C:\Python27\ArcGISx6410.4 is included? Is it a separate download and install? All I have been able to find is this link for 10.2.

In this link for 10.4, it says:

  Background Geoprocessing (64-bit) is available as a separate installation on top of ArcGIS for Desktop

But I can't find anywhere on that page that says how to go about the "separate installation"...


texts: Using Raster Bender plugin to rubbersheet raster layers in QGIS

texts:
I have a distorted raster of a historic map that I need to warp to match the modern map.
I've installed the plugin Raster Bender into QGIS 2.18, but I can't work out how to use it, and can find no instructions online.
How do I set up the pairs of points for it to use to warp the raster?
If there is an alternative method for QGIS 3, that would do. I have both versions.


texts: ArcGIS do calcuation on multiple input fields of attribute table and the number of input fields is changing

texts:
I have attribute table with headers:
ID(bgrp), F1 (F2011_46), F2 (F2011_48), ...Fn

I want to do math calculations on the fields from F1 to Fn. For example, count the total number of non-zero Fn fields for each ID unit. Normally, a field calculator using python script will do the job. But the problem is n vary each time and I have many such tables, so I don't know if I should modify the input fields each time when do field calculation for one such table. Additionally, n is usually large, 30 to 100, and I don't think it is efficient to manually select all the Fn fields in a Pre-Logic Script Code.


texts:
I have a huge raster file and want to increment all the values within a relatively small shapefile by 1 in Python. Speed is of critical importance, as I have to do this many times.

I have been looking at:
http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html

However, changing the code to my purposes, it seems elaborate to clip the raster, then set all values to one and update the raster with it's own clip.


texts: Drag and drop interaction error handling in OpenLayers

texts:
I'm trying to handle errors when doing drag and drop interaction for GPX files in OpenLayers 6.4.3, like for example when file being dropped has wrong format.
I found out that in the case of parse error there is now way to catch the error and inform user about that. The only visible consequence of error is this error in the console (but only for some browsers), which obviously comes from XML parser:

XML5619: Incorrect document syntax. Line: 1, Column 1

I tried to hook on interaction's error event, but it's not triggered at all.
Here is the relevant code:
var dragAndDropInteraction = new ol.interaction.DragAndDrop({
  formatConstructors: [
    ol.format.GPX
  ]
});

var map = new ol.Map({
  interactions: ol.interaction.defaults().extend([dragAndDropInteraction]),
  layers: [
    new ol.layer.Tile({
      source: new ol.source.OSM()
    })
  ],
  target: 'map',
  view: new ol.View({
    center: [0, 0],
    zoom: 2
  })
});

dragAndDropInteraction.on('addfeatures', function(event) {
  var vectorSource = new ol.source.Vector({
    features: event.features
  });
  map.addLayer(new ol.layer.Vector({
    source: vectorSource,
    style: styleFunction
  }));
  map.getView().fit(
    vectorSource.getExtent(), {
      constrainResolution: false,
      padding: [20, 50, 20, 50]
    });
});

Is there a way to catch parse error and consequently inform user about that?


texts: Using more than a single layer as barriers in arcpy.TraceGeometricNetwork_management()

texts:
I am working on script that iterates through trace flags and traces upstream.
For each trace flag, I have selected features from multiple layers in the geometric network that I want to act as barriers to terminate the trace. I have tried merging all of the selected features in each iteration to a single feature class, but if there is not one and only one feature at each junction (we have some points that are disconnected and coincident with connected junctions) then sometimes the disconnected feature is flagged as the barrier and the trace does not end.

I have instead made a list of layers that have selected features as barriers.

list_of_layers =['lyr1','lyr2']
arcpy.TraceGeometricNetwork_management(geomNetwork, r"outNet",traceflag ,"TRACE_UPSTREAM", list_of_layers, "", "", "", "", "NO_TRACE_ENDS", "NO_TRACE_INDETERMINATE_FLOW", "", "", "AS_IS", "", "", "", "AS_IS")

but this fails

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.7\arcpy\arcpy\management.py", line 6253, in TraceGeometricNetwork
    raise e
RuntimeError: Object: Error in executing tool

An earlier comment directed me to the tool help with states that this parameter is to be feature layer, not a list of barrier layers

How might one perform an upstream trace with barrier features selected on multiple point layer?


texts: Conversion from 3 dimension geometry to 2 dimension

texts:
I have a Point table in which the geometry column contains both 3 dimension and 2 dimension values.I want to convert this 3 dimensional geometry to 2 dimension.Is there any function in PostGIS? Please help me to solve this problem.

Regards
Sreesha.T.A


texts: How to aggregate pgr_dijkstra results by attribute?

texts:
I working with a table with a lot of points per road/street. I'm using pgRouting's pgr_dijkstra method to give me a path, but it returns too many points per road/street, I'd like to know if there is a built-in method in pgRouting to aggregate the results by column, in this case the same street name IF it exists. I want to transform the results table:

seq | cost |   street   |     st_asgeojson     |
----+------+------------+----------------------|
 1  |  2.5 |  Street 1  | {coordinates: {...}  |
 2  |  4.5 |  Street 1  | {coordinates: {...}  |
 3  |  2.1 |  Street 2  | {coordinates: {...}  |
 4  |  1.9 |  Street 2  | {coordinates: {...}  |
 5  |  2.9 |  Street 1  | {coordinates: {...}  |
 6  |  3.0 |            | {coordinates: {...}  |
 7  |  4.5 |            | {coordinates: {...}  |

into this:

seq | cost |   street   |     st_asgeojson     |
----+------+------------+----------------------|
 1  |  7.0 |  Street 1  | {coordinates: {...}  |
 2  |  4.0 |  Street 2  | {coordinates: {...}  |
 3  |  2.9 |  Street 1  | {coordinates: {...}  |
 4  |  3.0 |            | {coordinates: {...}  |
 5  |  4.5 |            | {coordinates: {...}  |

the current query is something like this:

SELECT seq, cost::double precision, street::character varying(192),
    ST_AsGeoJSON(geom)
    FROM pgr_dijkstra('SELECT gid AS id,
                              source::integer,
                              target::integer,
                              cost::double precision AS cost
                              FROM my_table,
                                  get_nearest_vertex_to_lon_lat(lat_st, long_st),
                                  get_nearest_vertex_to_lon_lat(-9.1511504, 38.7487786),
                                        false,
                                        false)
    a LEFT JOIN my_table b ON (a.id2 = b.gid);


texts: QGIS auto backup... is there one I can set?

texts:
I see QGIS makes a backup of the last manually saved project.  But how about a timed interval that automatically makes backups as I work.  Just had a QGIS crash dump, lost lot's of work.  Can't save after every item I edit.

If there is not one I am surprised this has not been requested as a feature in the past.


texts: Geoprocessing breaks my data links after upgrade?

texts:
I've had this problem all day, and it's happened before, too. I am running desktop 10.3.1 and want to produce a very simple buffer from a single parcel. The output is assigned to the same .gdb that the buffer output is assigned to. When I run the buffer, every other feature class in the map from that .gdb breaks the data link and I end up with a map full of red exclamation points, except for the most recent layer.
Occasionally this happens even when all I am doing is exported selected features into their own layer. Has anybody else experienced this?


texts:
NASA recently published this long imagery swath onto Gigapan.

Is there some way to mash up imagery from Gigapan so that it displays in a web map, using one of the ArcGIS web mapping API's?


texts: Earth Engine: Simplest way to move from EE image to array for use in sklearn?

texts:
I've produced a very large multiband image in EE with the goal of classifying it using the classifiers implemented in sklearn (the native ones implemented in EE don't provide enough flexibility for my purposes). sklearn uses 2-D arrays, so minimally I would need to convert each band to a 2D array and feed them in separately as explanatory variables. That's all fine.

Here's my problem: With a raster covering >150k km2, it is beyond tedious and cumbersome to Export.image.toDrive for each band, only to then re-import them to a python environment using rasterio. Ideally there would be some way to convert EE image objects to sklearn-readable NumPy arrays directly using the EE Python API (Google seems to tease as much with their documentation touting the advantages of using EE in Colab: "Seamless integration with Python data science libraries").

Is there a straightforward way to do this that I'm missing?


texts: Combining neighbor polygons (NOT touching) making one using JTS/GeoTools?

texts:
This question is similar to:

Combining touching polygons, keeping all attributes from Class "A" feature, using ArcGIS Desktop? and
Combining adjacent polygons sharing similar trait into single polygon using ArcGIS Desktop?

but need to do this programmatically with JTS.

This question is not for overlapping polygons or touching polygons but for adjacent polygons that is the difference when they don't touch but next to each other.
the answer to this question Remove intersect/touch from polygons, to making one with JTS/GeoTools , worked for intersecting..but not for adjacent.

I am trying to  implement this using JTS/GeoTools.. i found the .dissolve(..) thinking that would do the trick but it is not working.

The original GeoJson( //loadAdjacentGeoJson() )

https://gist.github.com/boundaries-io/ef7c423f3e174de0f185cb86a1e873e3

Updated Solution:

 for (GeoJsonPolygon geoJsonPolygon : geoJsonPolygonList) {
                Coordinate[] coordinates = toCoordinateArray(geoJsonPolygon);
                SimpleFeatureBuilder simpleFeatureBuilder = new SimpleFeatureBuilder(simpleFeatureType);
                simpleFeatureBuilder.add(inflate(geometryFactory.createPolygon(coordinates)));  //add inflated polygons

                SimpleFeature simpleFeature = simpleFeatureBuilder.buildFeature(Integer.toString(index));
                simpleFeature.getUserData().put("zipcode", postal.getZipCode());
                index = index + 1;
                geometryList.add((Geometry)simpleFeature.getDefaultGeometry());
            }

        Object obj = geometryFactory.buildGeometry(geometryList);

        GeometryCollection geometryCollection = (GeometryCollection) obj;

        Geometry union = geometryCollection.union();  //perform union

        LOGGER.info(new GeoJsonWriter().write(geometryCollection)); //shows inflated polygons with overlap

        LOGGER.info(new GeoJsonWriter().write( deflate( union ))); //deflate back.

inflated polygons so union will work.:

https://gist.github.com/boundaries-io/80614499880bb4dabb57205cc529ca7c

inlate/deflate methods:

private Geometry deflate(Geometry geom) {
        BufferParameters bufferParameters = new BufferParameters();
        bufferParameters.setEndCapStyle(BufferParameters.CAP_ROUND);
        bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
        Geometry buffered = BufferOp.bufferOp(geom,-.0001, bufferParameters);
        buffered.setUserData(geom.getUserData());
        return buffered;
    }

    private Geometry inflate(Geometry geom) {
        BufferParameters bufferParameters = new BufferParameters();
        bufferParameters.setEndCapStyle(BufferParameters.CAP_ROUND);
        bufferParameters.setJoinStyle(BufferParameters.JOIN_MITRE);
        Geometry buffered = BufferOp.bufferOp(geom,.0001, bufferParameters);
        buffered.setUserData(geom.getUserData());
        return buffered;
    }

results:


texts: Postgis 2: Polygon/MultiPolygon in a single column

texts:
Storing MultiPolygon and Polygon features in a single column used to be as simple as not creating or removing the geometry type check constraint.

Is there a similar simple fix in PostGIS 2?


texts: How do I add a space in a label with vbscript?

texts: I am labeling point features in Acrmap 10.3.  I have the points labeled with three fields from the attributes on the same line, but there is no space between them.  what is the expression to add a space?


texts:
I'm trying to reorder the TABLES in my ArcMap layer list. I can't find a way to do it. I can reorder normal layers since they appear in the "List by Drawing Order" tab of the "Table of Contents".

The reason why the order is important for the tables is that when I publish this to an ArcGIS Server as a map service the tables get an ID based on their order. If their ID changes in any way then code which looks them up based on that ID will break.

The question is, in ArcMap, can I reorder my TABLES in my layer list at all?

ArcMap 10.2
Oracle 10g


texts: Georeference a dataset in Carto with imported data

texts: I have a number of datasets (in CSV format) that I would like to upload into Carto - the issue is that the dataset with the postal code/georeference data is separate. I would like to upload my data and then use the dataset with georeferenced data (postal code) to display on a carto map.


texts: How to change dialog box fonts (a tool created with ModelBuilder)

texts:
I have created a model using ModelBuilder and now is the time to publish the model. But the problem is that I need to change the fonts of the dialog box. Does anybody knows how can I change these?
P.S. I have tried to change the Stylesheets, but nothing has changed.


texts: Function with map.addLayer() and map.removeLayer() not working

texts:
I am trying to create buttons that add and remove layer groups and I am not having any luck. Can you take a look at how I have this set up?
    var main = L.tileLayer('https://api.mapbox.com/styles/v1/{id}/tiles/{z}/{x}/{y}?access_token={accessToken}', {
            attribution: 'OpenStreetMap, CC-BY-SA, Mapbox',
                maxZoom: 18,
                id: 'mapbox/streets-v11',
                tileSize: 512,
                zoomOffset: -1,
                accessToken: 'pk.eyJ1IjoibWljaGFlbGNwIiwiYSI6ImNrN3V5eHVwNTAxMzgza2x1d2s1N3lkbmIifQ.Rmgw9tQHXVTL5_w6jvGjlw'
        });

    var map = L.map('mapid', {
                 center: [59.3307804,18.0569225],
                 zoom: 11,
                 layers: [main],
                 fullscreenControl: true,
                 zoomControl: false,
                 dragging: !L.Browser.mobile
        });

var zoomHome = new L.Control.zoomHome();
zoomHome.addTo(map);

var baseMaps = {
       "All Activities": allWater,
       "Beaches": allBeaches,
       "Wading pools": waterR,
       "Swimming pools": pools
  };

var control=  L.control.layers(baseMaps, null, {collapsed:false});
control.addTo(map);

Here are the functions I have written but they don't seem to work. I am using onclick= to trigger the functions
    function WaterFilter0(){
  map.addLayer(allWater)
    if(map.hasLayer(allBeaches)){
        map.removeLayer(allBeaches);
      } else if(map.hasLayer(waterR)){
        map.removeLayer(waterR);
      } else if(map.hasLayer(pools)){
        map.removeLayer(pools);
    }
}

function WaterFilter1(){
  map.addLayer(allBeaches)
    if(map.hasLayer(allWater)){
        map.removeLayer(allWater);
      } else if(map.hasLayer(waterR)){
        map.removeLayer(waterR);
      } else if(map.hasLayer(pools)){
        map.removeLayer(pools);
    }
}

function WaterFilter2(){
  map.addLayer(waterR)
    if(map.hasLayer(allWater)){
        map.removeLayer(allWater);
      } else if(map.hasLayer(allBeaches)){
        map.removeLayer(allBeaches);
      } else if(map.hasLayer(pools)){
        map.removeLayer(pools);
    }
}

function WaterFilter3(){
  map.addLayer(pools)
    if(map.hasLayer(allWater)){
        map.removeLayer(allWater);
      } else if(map.hasLayer(allBeaches)){
        map.removeLayer(allBeaches);
      } else if(map.hasLayer(waterR)){
        map.removeLayer(waterR);
    }
}
</script>

This is code is live here: https://www.leksplay.com/tabs


texts: How to use GRASS GIS in Python standalone?

texts:
I'm trying to use GRASS module in a Python Skript. I already checked this:
Is it possible to use GRASS GIS in Python stand alone scripts?
and this Using Python script to control GRASS GIS from outside?
but it doesn't help with my problem.
I did pip install grass-session which worked well. My script is just:
import grass_session and when I run it I get the error RuntimeError: Cannot find GRASS GIS start script: 'grass', set the right one using the GRASSBIN environm. variable.
First I tried GRASSBIN="$path\\to\\software\\grass\\grass78 python grass.py" as suggested in
https://grasswiki.osgeo.org/wiki/Working_with_GRASS_without_starting_it_explicitly in the block Python: GRASS GIS 7 with an external library: grass-session but I still got the same error.
After checking http://osgeo-org.1560.x6.nabble.com/GRASSBIN-environment-variable-td5413301.html I solved this error adding the following line before the import line: os.environ['GRASSBIN']=r"path\to\software\grass\grass78", although I feel this only should be necessary if my grass version is different from 7.8.
Now when I run these two lines first I'm happy to see that I don't get the error again but after a while I notice that nothing happens. It keeps running without any exit code unless I stop the script manually. I would expect the script to import the module and finish with exit code 0.
I have grass gis 7.8.2. installed with QGIS and Python 3.7 and working on Windows 10 Pro 64 bit.


texts: Cause of selected [ArcSDE 10] feature dataset does not contain any feature classes which can participate in topology?

texts:
When trying to create a topology for a polygon feature class in ArcSDE 10.  I right click the feature dataset that the feature class resides and try to create a new topology.  I get an error saying, "The selected feature dataset does not contain any feature classes which can participate in a topology".

The feature class I am targeting has a z value.

Would this be the culprit?


texts:
I am having troubles with extracting data in ArcGIS buffer extractions from the MOD11A2files after processing them into GEOtiff, most of them seem fine but when looking at nightly temperature on some of my coordinates of interest my values go way below zero. Some down to -273. I have multiplied the data with 0,02 and then subtracted 273 to get them as Celsius. Is this a common problem? I can't manage to find anything about it in forums etc and unfortunately I have a tight schedule and am stressed to get the results out.

Anyone familiar with this problem?


texts:
I just tried to use Reports in ArcGIS Pro 1.2, and have not been able to find them in either the application or its documentation.

A GeoNet thread asked the same question for ArcGIS Pro 1.1 but that has not been replied to.

Has anyone found them in ArcGIS Pro 1.2, or know when they are expected to be added?


texts: Fastest way to write large numpy arrays to feature class: as rows, or as columns?

texts:
I have a numpy data processing loop that repeats ~50,000x and outputs ~23 data values. I want to write this data to the attribute table of my ~50,000-point multipoint feature class.
With each loop, I generate 8 individual values, plus 5 sets of three values. (I'm still deciding whether to group them as 3 sets of 5 or 5 sets of 3, TBD by the answer to below.)

With this number of iterations, I want to be mindful of processing time. How should I approach organizing & writing the data? Should I:

...write 23 fields of data to the output file row-by-row at the end of each loop? (Or is the processing time of writing to file too long to do ~50Kx?)
...or collect them into 23 separate ~50K-length numpy arrays and write them column-by-column after all looping is done? (Or would I likely run into memory issues?)

Here's my approach for the second method:

#23 data columns to be added
fieldlist = ["azimuth", "orient", "crit#f", "freq1", "amp1", "pwr%1", "freq2", "amp2", "pwr%2", "freq3", "amp3", "pwr%3", "freq4", "amp4", "pwr%4", "freq5", "amp5", "pwr%5", "Nspurs", "len_mean", "len_sd", "ht_mean", "ht_sd"]

for fieldname in fieldlist:
    arcpy.AddField_management(gridcopy, fieldname, "FLOAT") #Add new fields
    with arcpy.da.UpdateCursor(gridcopy, [fieldname]) as cursor: #Add data
        i = 0
        for row in cursor:
            row = X[i]
            cursor.updateRow(row)
            i += 1


texts: Understanding geography, geometry, distance and wkt in PostGIS

texts:
I'm trying to understand the data conversions that happen in PostGIS when calculating distances in SRID 4326.

Why is it that when I run on my marin_jobs points table:

select cast (the_geom as geography), the_geom from
marin_jobs where blockid10 = '060411141001010' limit 1;

I get

the_geom      the_geom
geography     geometry
010100....    010100....

And the wkt output values are exactly the same?

I came to this question because when I run a distance query (see below for code) on two points in my dataset, if I don't use cast () as geography I get out degrees, just like the Boundless Geography tutorial said I would. http://postgis.net/workshops/postgis-intro/geography.html

But I don't understand how if the inputs are (seem to be anyway) exactly the same, the outputs are different. Or should I ask...are the outputs really different?

Does that mean that the value of 0.0616262566763201 degrees and 6007.875564598 meters are the same, but just in different units?

But this doesn't jive with what the Boundless Tutorial says,

"On a sphere, the size of one “degree square” is quite variable, becoming smaller as you move away from the equator. Think of the meridians (vertical lines) on the globe getting closer to each other as you go towards the poles. So, a distance of 121 degrees doesn’t mean anything. It is a nonsense number."

Any plain English claficiation would be really helpful.  This is keeping me up at night.

Distance Queries:

select st_distance (
               (select the_geom from marin_jobs where blockid10 = '060411141001010' limit 1),
               (select the_geom from marin_jobs where blockid10 = '060411070005007' limit 1)
               );

output: 0.0616262566763201

select st_distance (
                (select cast (the_geom as geography) from marin_jobs where blockid10 = '060411141001010' limit 1),
                (select cast (the_geom as geography) from marin_jobs where blockid10 = '060411070005007' limit 1)
               );

output: 6007.875564598


texts: Customize behaviour of Ol.interaction.draw.createbox() geometryFunction?

texts:
I'd like to enable a draw function in an Openlayers 3 application for drawing a rectangle - up to a certain size (in the defined projected coordinate system, in square meters).
For this ol.interaction.draw instance the draw.createBox geometryFunction would be used and the 'Circle' type, as per

http://openlayers.org/en/latest/examples/draw-shapes.html?q=draw

the geometryFunction is where to "modify the geometry during drawing" as per a prior thread:

Undoing last point when drawing linestring in OpenLayers 3? )

but for drawing a rectangle the predefined createBox geometryFunction is used and not a custom one that might perhaps continually check the size of the rectangle. Or is this not the place to put in a check of the current size of the rectangle being drawn? Is what I wish realistic?


texts: How to connect QGIS to an online database on Mac OS X?

texts:
I'm a real beginner/dummy when it comes to QGIS and online databases. I have tried a variety of things trying to connect to Spatialite, PostGIS and MySQL but all attempts have failed due to my lack of knowledge on the subject.

So I was wandering:

What is the easiest way to connect QGIS to an online database?
Is there such a thing as a step by step manual to do this that is easy to understand by non programmers and dummies like myself?

Many thanks for your help

Ed

specs:
MAC OS X Mountain Lion (but I also have Ubuntu)

QGIS 1.8 Lisboa


texts: How to ascertain if a cell/field contains a value in ArcGIS attribute table

texts:
I have done a bit of research on this but can't seem to find an answer. I am wondering if there is a way to identify whether a cell/field in the ArcGIS attribute table contains a value using the search cursor. Initially, I thought it would be as simple as this:

arcpy.MakeFeatureLayer_management(node110, "node_layer12")
arcpy.SelectLayerByAttribute_management ("node_layer12", "NEW_SELECTION", query1)

rows = arcpy.SearchCursor("node_layer12")
for row in rows:
        link_id2 = row.getValue("DRN_link")

if link_id2 == "E*":
    print "YES"

unfortunately this didn't work. I knew that if the field I was searching had a value it would start with an 'E' so I thought a wildcard might work. Obviously not. Any ideas?

Please note. I only have access to arcgis 10 so can't use the arcpy.da.searchcursor.


texts: OpenLayers Plugin Error: latitude or longitude exceeded limits

texts:
I'm aware the 'traceback' OpenLayers plugin error has been dealt with before, so apologies for raising the issue once again. Unfortunately for me, the solutions posted were way over my head, e.g. using workarounds and code posted on github, etc. I find it all mind-boggling.

If anyone could post a step by step answer for an absolute beginner like me I would be eternally grateful! Thank you. I hope that one day I can contribute something.

The error message reads:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/MCMA1580/.qgis//python/plugins\openlayers_plugin\openlayers_plugin.py", line 48, in addLayer
    self.__plugin.addLayer(self)
  File "C:/Users/MCMA1580/.qgis//python/plugins\openlayers_plugin\openlayers_plugin.py", line 190, in addLayer
    self.__setMapSrsGoogle()
  File "C:/Users/MCMA1580/.qgis//python/plugins\openlayers_plugin\openlayers_plugin.py", line 249, in __setMapSrsGoogle
    extMap = coodTrans.transform(extMap, QgsCoordinateTransform.ForwardTransform)
QgsCsException: forward transform of
(0.884883, 2.549197)
PROJ.4: +proj=longlat +datum=WGS84 +no_defs +to +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs
Error: latitude or longitude exceeded limits


texts: Getting maximum value of measured values every xx meters

texts:
I have huge CSV files from measurements (most of them are done by car), like this:
measurement file
The measurement system record one value every second and save this together with the GPS coordinates. That's why there are a lot of measure points very close together (e.g. because of waiting at traffic lights or a part of the streets are run and measured multiple times. In the end, the file is bigger than necessary.
Now I'm looking for a solution to reduce the measurement points.
At the end of the day, I want to have a measured value displayed on the map only every 5 meters if the measured positions were closer than 5 meters. If the measured values were further apart, then of course these are the way they are.
If the points are closer together than 5 meters, the new center point coordinate should contain the values of the measured value that contains the maximum value of the RSRP column within these 5 meters.
Most important is for me, to have the maximum value of RSRP in a 5m radius at the end.
If the assigning of a new 5m center coordinate cannot be easily done with QGIS tools, it would also be sufficient for me if the points with the maximum value of RSRP are displayed within 5m at their original position = the points do not have a fixed 5m distance, but between 0m and 10m distance.
While the simplification, the points should not be summarized with their chronologically adjacent points only, but with all in 5m radius (those that were recorded at a later point in time by repeatedly measuring the same route section)
With a solution from @Taras, I could remove a lot of intermediate positions and make the file smaller. However, this solution has the following disadvantages for me and I am looking for a solution that also solves this:

With 'Extract by location', the point value of the base point of the simplified path is adopted (regardless of the RSRP value). But I absolutely need the maximum value of RSRP in a 5m radius.
The resulting points are at the end not with a fixed 5m distance. This means that if the path is straight for a long time before the simplification, the path will not have any support points over long distances after the simplification to which the original measuring points would be mapped in the next step. So I lose the information which measured values the points (5m apart from each other) have on the way in between. If this cannot be easily solved with QGIS tools, it would also be sufficient for me if the points with the maximum value of RSRP are displayed within 5m at their original position = the points at the end do not have a fixed 5m distance, but between 0m and 10m distance.
In the simplification, the points are only summarized with their chronologically adjacent points, but not also with those that were recorded at a later point in time by repeatedly measuring the same route section (cf. in Rathen the route between the ferry and the intersection)


texts: I have applied a point clustering effect on a layer in QGIS. What I'd like to be able to do is cluster on all features except the ones that I have a NULL value in an attribute. I can Categorize on feature when zoomed in using the rendered Settings, but there doesn't seem to be an option for categorizing data when clustered.


texts: Join attributes from point into lines (but only lines with certain attribute value)

texts:
I have point objects with attributes and line objects, where the attributes are missing.
The only connection between these geometries is an angle I calculated. (the labeling of the points has the same angle than the lines)
So I want to join the attributes from the points Into the line object, but only in the line object with the same angle value as the point. (There is a lot of overlap because I have over 10,000 line objects and each point belongs to a specific line)
My idea is it to generate a buffer around the points which hits the lines. Next step would be to check if the angle value matches before joining the attributes and this part is where I'm struggling.
Is there a possibility to do this?
I'm new to GIS.


texts:
I've seen many ask how to clip a raster via a polygon mask but have yet to find the opposite operation. How does one clip a polygon layer via raster?
ta
Specifically how does one clip a polygon layer via certain raster values?

For example imagine a flooding model. I would like to get the polygon shape that represents a "safe" area in a certain region. I would have to clip the raster extent against a certain elevatio value in the raster.

Any ideas?


texts:
Is it possible to, through Python, to create the relations that it's possible to create in Layers/Properties/Fields/Relations, and if so how?
I first thought these were saved in the styles, but that was not the case.


texts: How to limit viewing to map boundary on Android?

texts:
I have an offline map of a city displaying on Android. Since my map is just small, when you scroll beyond the the boundaries of the map you see only a gray screen and you will probably be lost and could not be able to view the city map.

How can I limit the viewing of the map to the boundaries of the city map? What should I do for the user to only scroll on the map area, and not beyond it?

My resources are OSM Mapnik, osmdroid-android-3.0.1.jar library.

Thank you for your time.


texts: Implement r.cross / ArcGIS Combine with GDAL

texts:
I am using GDAL to process a lot of rasters. I then need to combine all the arrays to create a single output raster that represents all the unique combinations of values.
I.E the same thing as r.cross in GRASS or Combine_sa in arcgis.

I'd like to program a similar thing using the GDAL library.
What is the algorithm? How would I handle rasters that are too large to fit in RAM? I.e., when chunking them into smaller arrays, how do I ensure that a combination in the nth chunk is given the same value if it appeared in a different chunk?

I'm struggling to come up with anything that will perform well.


texts: How to implement Google Earth style zoom in/out and pan animation in ArcGIS Engine applications?

texts:
I am developing an application using ArcGIS Engine 10 and the C# language.

I wish to add some more animation capabilities into my application.

In other words, I want to add animation to zoom in, zoom out and pan like the animation in Google Earth's pan and zoom in/out.

I use OpenStreetMap (OSM) maps for the base map and also have other local layers in my application.

Is there any way to use animations like Google Earth animation for zoom in/out and pan to act with my map ?


texts: obtaining the viewer extent of an ArcGIS portal page w/ python widget

texts:
I will be hosting various locations of features on an ArcGIS portal page, and want to add a geoprocessing widget (Python toolbox) that will export those features currently seen in the viewer.

The widget will select features that exist within the current view, so I want to grab the coordinates of the viewer extent. I've searched extensively for some ArcPy/Python code that will simply get the current extent of the ArcGIS Portal viewer at the time the user clicks the widget... So they do not have to set a polygon or circle to grab features, it will simply record the extent coordinates of the viewer.

Any thoughts on how this can be accomplished? Is there an ArcPy command that is escaping me?


texts: Heat Map- how to add colour division and legend

texts:
I would like to ask if anyone knows how to colour the certain counties of England + legend about the colour's meaning. Here is my code which will give me the chosen counties of England on a map. But I would like to add numbers from my excel file about how many repairs on which properties are on each county similar to the link below.

GB0 <-readRDS("gadm36_GBR_0_sp.rds")
GB1 <-readRDS("gadm36_GBR_1_sp.rds") # England, Ireland, Scotland, Wales
GB2 <-readRDS("gadm36_GBR_2_sp (1).rds")
GB3 <-readRDS("gadm36_GBR_3_sp.rds") # Cities

x<-(GB3$NAME_2=="Wirral" | GB3$NAME_2=="Sefton"|GB3$NAME_2=="Lancashire"|GB3$NAME_2=="Cheshire East"|
      GB3$NAME_2=="Blackpool"|GB3$NAME_2=="Blackburn with Darwen"|GB3$NAME_2=="Rochdale"|GB3$NAME_2=="Bolton"|
      GB3$NAME_2=="Oldham"|GB3$NAME_2=="Wigan"|GB3$NAME_2=="Salford"|GB3$NAME_2=="Knowsley"|GB3$NAME_2=="Warrington"|
      GB3$NAME_2=="Liverpool"|GB3$NAME_2=="Manchester"|GB3$NAME_2=="Stockport"|GB3$NAME_2=="Trafford"|
      GB3$NAME_2=="Halton"|GB3$NAME_2=="Tameside"|GB3$NAME_2=="Cheshire West and Chester"|GB3$NAME_2=="Bury"|
      GB3$NAME_2=="Saint Helens")

County<-GB3[x,]
plot(County)
plot(County[x,],border="blue",add=TRUE) #Notice add=T; this draws on top of current plo
text(coordinates(County),
     labels=County$NAME_3,col="red",cex=0.7)


texts: How to hide grid's labels when zooming out?

texts:
I have used polygon grid layer on my map, and I used to label the layer with ID automatically, and the result was convenient for me:

But when I zoomed out, the labels didn't hide at all! and they keep appearing in frustrating way!

I would like to get the result when I zooming out like this:

What to do in this case to get this result?
How to hide grid's labels when I zooming out?

My vector grid is 1000 meter polygon.


texts: Browse (add/remove) ArcGIS Server Services (Layers) from the ESRI Javascript API viewer

texts: I am quite new to the JavaScript API. I am not if there is a widget that does what I would like but I am trying to create a web map that allows the users to choose/add/remove layers located on our GIS Server. I can't figure out how to achieve something similar as the picture attachment. If anyone have any samples/tips, please let me know.


texts: Correctly display the colors in MATLAB colormap and colorbar?

texts:
I am plotting spatial maps with fill values ranging from 0.05 to -0.05

Below is the part of code for plotting colormap and colorbar

    colormap(brewermap([10],'RdYlBu'));
    caxis([-0.05 0.05])

and here is the figure

Here the problem is, each map has majority of the 0 fill values, but the color map show them in between 0 to 0.01 as depicted through colorbar.

Is it possible to assign some other color to 0 values?


texts:
Some polygons in my map have a "name" (but no "capacity") and some have a "capacity" (but no "name") and some have both.

I want to label the polygons with the "name" and the "capacity". Therefore I use the following string:

"name" || '\n' || "capacity"

But this does not work.

Can someone help me?


texts: Problem by changing the order of Layers in QGis Layer Tree

texts:
I'm using QGis 2.4 on a virtual machine (Ubuntu 14.4 with XFCE). I'm connected to the VM via vncviewer, what works pretty well.

But now I want to add some more layers to an existing QGis Project. After doing that, I can't change the order of the layers by Drag And Drop!

Is there anybody who knows that? Any Idea how to solve this problem!


texts: Save layers from a gdb as sf objects in a loop

texts:
I'm taking an R class and I have an assignment that I'm trying to figure out. Instructor wants us to save all layers in a gdb file as sf objects. He told us that we should use loops.
I'm able to import the gdb file using
scd <- st_read(dsn='/path/to.gdb')

And I'm also able to list layers using
scdl <- st_layers(dsn='/path/to.gdb')

Update:

  st_read(dsn = '/path/to.gdb', layer=i)
  feat <- st_read(code here)
}

i
this is what I have done so far which loops over the layernames. but when I try to assign them to objects using feat <- st_read(code here) line, it breaks my loop.
"Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet,  :
Expecting a single value: [extent=5].
In addition: Warning message:
In if (nchar(dsn) < 1) { :
the condition has length > 1 and only the first element will be used"


texts: Can you download all google maps data in certain area?

texts:
Knowing the reference point, I want to get a list of all the shops (and co ordinates) based in that area.. For example all the shops in this mall:
https://goo.gl/maps/GFfWx8pze2T2

Is there a tool to do this online or something I can programatically use?


texts: Using folder wildcards to copy files with Python

texts:
I would like to copy the data of one specific shapefile from a specific folder ("Grid") to a number of folders containing the name "Country". The country shapefiles would keep their original name. Since there are multiple Country folders, I thought using a wildcard would be quite useful. The following diagram tries to illustrate what I am doing:

Below is the code I have at the moment where I used this post as a guide:

import os, sys, glob, shutil

root_dir = "C:\Users\xxxx\Desktop\Test\\"
country_dir = "Country*\\"

#do_some_function

for path, dirname in os.walk(root_dir):
    if country_dir in path:
        for subdirname in dirname:
            if subdirname.startswith('Country'):
                for fname in glob.glob('*.shp'):
                    shutil.copyfile(do_some_function, fname)

However, I receive the following error:

  ValueError: too many values to unpack

Any advice on what I am doing wrong?

EDIT:

Following @NathanW's suggestion, I've changed the script to the following:

import os, sys, shutil

root_dir = "C:\Users\xxxx\Desktop\Test\\"
country_dir = "Country*\\"

#do_some_function

for path, dirs, files in os.walk(root_dir + folder_dir):
    for file in files:
        if file.endswith('*.shp'):
            shutil.copy(do_some_function, file)

Unfortunately, nothing happens.


texts: arcpy.da.Editor not saving edits to base Oracle Spatial tables

texts:
Can't post code because of Fed policy, so I will ask the question here.

We are working on a series of polygon feature classes that are used to define large activity areas. These polygons are modified using arcpy such that activity exclusion areas are cut out of the larger polygon. To accomplish this end we are using arcpy.da.Editor(). The code works beautifully in a file geodatabase, but when executed in Oracle Spatial, it refuses to push edits from the delta tables into the 'A' table in Oracle Spatial. My first thought to a solution here would be to access the 'SHAPE@' object in the fGDB then convert it to a string and use an "UPDATE" command in Oracle SQL via ArcSDESQLExecute to push the edits to the 'A' tables.

Can anyone else come up with a better solution? Our DBA insists that this issue is on the code side and not a quirk of Oracle.


texts:
I georeferenced some things about a month ago no problem. However now when I click my reference plugin it does not open the Georeferencer dialogue box or the georeference drop down menu. It is enabled it just does not start.

Has this happened to anyone?


texts: How to derive a hydrologically correct DEM from point shapefile using QGIS?

texts: I used grass gis v.surf.rst function to interpolate a surface using vector points. But it is very difficult to find a optimum parameters to be used. I am still unable to develop a accurate terrain


texts:
I'm trying to use/learn some PHP MapScript on a GDB file and cannot get anything from it.

The PHP below works great if I use countries.sqllite (via countries.map) file I have and shows me an array of countries both in a for and while loop.

However, when I switch and point to a buildings.gdb file (via buildings.map) I get an error of "msQueryByRect(): No matching record(s) found." for the first for loop and blank or no data for the while loop.

Both .map files just point to the one layer and show fine as WMS in OpenLayers and in QGIS.

I've been searching through PHP MapScript API stuff all afternoon (starting with https://mapserver.org/mapscript/php/phpmapscript.html) and found nothing to help me, and am asking if there is anything I missed or any special considerations I need to use MapScript a GDB file.

$map_path = "/ms4w/maps/";
$map = ms_newMapObj($map_path."buildings.map");
$map = ms_newMapObj($map_path."countries.map");
print_r($map);

$layers = $map->getAllLayerNames();           print_r($layers);
$myLayerName = $layers[0];                    echo '|'.$myLayerName.'|
';

$l = $map->getLayerByName($myLayerName);      print_r($map->extent);
$l->queryByRect($map->extent);
for ($i=0; $i<$l->getNumResults();$i++){
    $s = $l->getShape($l->getResult($i));
    echo $i.'=';                              print_r($s->values);
}

$layer = $map->getLayer(0);                 print_r($layer);
$status = $layer->open();                   print_r($status.'
');
$extent = $map->extent;                     print_r($extent);
$status = $layer->whichShapes($extent);     print_r('
'.$status.'
');
while ($shape = $layer->nextShape()) {
    echo $shape->index ."=>";               print_r($shape->values);
}
$layer->close();


texts: Using Copy (Management) from GetParametersAsText() in ArcPy?

texts:
Having some issues getting the arcpy.Copy_management to work.
Essentially I simply would like to be able to have a user select a folder which contains multiple shapefiles and copy them to another user selected fodler without hard coding the locations and filenames. Both folders would be obtained by the GetParametersAsText() function.

At this stage all I want to do is get the copy code correct and make a copy of each of the datasets in the dataset folder. I may need to check for differing things such as spatialReference but for now getting the copy code sorted would be great!

I have set my workspace to the dataset directory using:

arcpy.env.workspace = datasetFolder

datasetFolder = arcpy.GetParameterAsText(0)
outputFolder = arcpy.GetParameterAsText(1)

The datasetFolder contains 3 shapefiles and I'm a bit lost on the parameters to use to copy each of them to the outputFolder directory. It surely isn't as simple as:

arcpy.Copy_management(datasetFolder, outputFolder)

All the files have been put into a list called inputFC.

Should the above Copy_management code work or is there a particular way to code the copy line so each files gets copied over?


texts: "buffer" created with DEM analysis - how to avoid it?

texts:
I have, to start, a DEM raster and a shapefile, with random created points, within a given polygon. The DEM raster was clipped with that same polygon.

With the DEM raster, I used GRASS to perform neighborhood analysis (using "r.neighbors"), and I applied a buffer of 150m. Since the pixel size is 30x30m, the size I used was 5 (and I did not select the "circular neighborhood" option).

Then, I uploaded the new raster to QGIS, to perform DEM Terrain models analyses (Raster > Analysis > DEM (Terrain Models)). In "Mode" option, I used: Slope, TRI and Roughness.

My final objective is to atribute the value of "Slope", "TRI" and "Roughness" to each point of the shapefile.
However, after creating "Slope", "TRI" and "Roughness", a kind of buffer (if it can be called this way) appeared surrounding the limits of the area, as you can see here:

The problem is, that buffer has an extremely high value that is not representative of the area, and some points fall into that buffer, so I obtain erroneous "Slope", "TRI" and "Roughness" values.

Do you have any idea of a solution for this problem?


texts: How to retrieve coded values (domain) from Runtime geodatabase?

texts:
I am working with the ArcGIS Runtime SDK for iOS 10.2.5 (not the latest version) in order to create an iOS app, using swift coding.

I have created a runtime geodatabase which includes coded values and I cannot find a way to retrieve them.

I checked the ESRI documentation for this version but I could not find something about it. Any ideas?


texts: Why is the CSV Layer converting a Date-Time field to Unix Timestamps?

texts:
I am using the ArcGIS JavaScript API to add a dynamic map in my application. I am loading data through a csvLayer.

The data looks something like this:

Here's a snippet of me initiating the csvLayer:

var csvLayer = new CSVLayer({
                    url: numberLayers[0].Url,
                    popupTemplate: popupTemplate,
                    featureReduction: {
                        type: "selection"
                    },
                });

The Popup template where I am displaying the information from the CSV file is given here:

var popupTemplate = {
                    title: "{SensorTitle}",
                    content:
                        [
                            {
                                type: "fields",
                                fieldInfos: [
                                    { fieldName: "DateTime", label: "Date" },
                                    { fieldName: "Temp", label: "Temperature (°F)" },
                                    { fieldName: "Humidity", label: "Humidity (%)" },
                                    { fieldName: "DewPoint", label: "Dew Point (°F)" },
                                    { fieldName: "PM1", label: "PM 1.0 (µg/m3)" },
                                    { fieldName: "PM25", label: "PM 2.5 (µg/m3)" },
                                    { fieldName: "PM10", label: "PM 10.0 (µg/m3)" },
                                    { fieldName: "O3", label: "Ozone (ppb)" },
                                    { fieldName: "NO2", label: "Nitrogen Dioxide (ppb)" },
                                    { fieldName: "CO", label: "Carbon Monoxide (ppb)" },
                                ]
                            },
                        ],
                    actions:
                        [
                            { title: "More Information", id: "more-info", className: "esri-icon-notice-round" }
                        ],
                    overwriteActions: true,
                };

Now this is how the data is presented on my application:

If you notice the datetime field for some reason changes automatically from a normal timestamp format to Unix epoch times. Why is this happening? How can I fix this?

P.S: If this post is not clear enough, just compare the two images I posted. The datetime field is automatically changing into a different format.


texts: How to get default QGIS context menu back after customizing it once?

texts:
I have created custom legend Layer context menu using

    class MyMenuProvider(QgsLayerTreeViewMenuProvider):
        def __init__(self, view):
            QgsLayerTreeViewMenuProvider.__init__(self)
            self.view = view

        def createContextMenu(self):
            if not self.view.currentLayer():
                return None
            m = QMenu()
            if type(self.view.currentLayer()) == QgsVectorLayer:
                m.addAction("Show Attribute Table",self.openAttributeTable)
            elif type(self.view.currentLayer()) == QgsRasterLayer:
                pass
            return m

       def openAttributeTable(self):
           self.attributeTable=attributeTableDialog()
           self.attributeTable.show()

but want to get the default context menu back (already provided by QGIS)!


texts: Convert spatial Geographic data from sql server 2008 to Json

texts: Is there any example to convert spatial data (point,line,polygon) stored as geographic in sql server 2008 into Json and display data using Arcgis Javascript Api 3.0.I have tried converting geographic data into Json given Geometry to Json  but as my data stored in  geographic data  format it failed.Thank you in advance.


texts: Using ArcPy GetCount to return number of records for each table in GDB or MDB

texts:
This is sample code from Esri
import arcpy
arcpy.env.workspace = "C:/data/data.gdb"
arcpy.GetCount_management("roads")

My understanding is this provides the number of records in the table for just the "roads" feature class.
I'm hoping to return the number of records for each feature class and table in the workspace.
Is there a simple way to do so, preferably in ArcCatalog?


texts:
Could anyone tell me how to select a date in the Symbology > Categories > Value field to express symbology for any cells that contain a date rather than having to add every date as a field to define my symbology?

I'm running this data from an excel sheet that gets continually updated so can not pick individual dates, I'd just like to know if there's a certain word or code for "any date" or "any value" that will trigger the symbology to work?


texts: Does gdalwarp resample when calculating geoid correction?

texts:
I wanted to perform a geoid offset correction to a DEM and saw this really useful post: Vertical Datum conversion of a raster with gdal

I noticed that although they use gdalwarp, they do not specify a resampling method (-r). Does gdalwarp always resample the image, making it important to specify the resampling method, or does it not resample when performing geoid correction? As the default resampling method is NN, I would like to avoid this when using elevation data but I don't want to force a resampling method if one isn't required.


texts: Self organizing maps and spatially Constrained Clustering in R

texts:
There are a few algorithms to perform spatial clustering in R but most of them focus on point pattern analysis (for example here). I would like to perform spatial constrained clustering for polygon data, something similar to this. There are some interesting Open-Source Tools to do this task like SOMVIS and other tools on the site from this link.
But rather than using "another" tool I would like to try this in R since I work a lot in R and I'd like to keep that flexible working environment.

So there seems to be very good information here. Does anybody have expirience with a similar task and can anybody give me some advice on good ressources, packages, books etc?

Best,Johannes


texts: GDAL or rasterio to find minimum extent of multple files in command line

texts:
I have two files, GeoTIFFs (.tif) in this case which need to be cropped to each other's minimum overlap, and have matching spatial resolutions. I also want to round their Pixel Size to a whole number (50m in my case). I strongly prefer using gdal or a similiar alternative like rasterio.

File1:

Upper left (240650.486, 4234958.562)
Pixel Size (50.000101..., -50.000101...)

File2:

Lower Left (242715.000, 4234848.000)
Pixel Size (50.00023..., -50.00023...)

Say I found the min overlapping extent in python for instance:

left_bounds = max(min(xcoords1), min(xcoords2))
right_bounds = min(max(xcoords1), max(xcoords2))
#same for lower_bounds and upper_bounds

I still couldn't simply gdalwarp -te left_bounds lower_bounds right_bounds top_bounds because my resolution still needs rounding.

Essentially I would like to do all this through gdal commands, but am curious the most parsimonious string of gdal commands that could clip to bounds, ensure alignment of two tifs, and round resolution to 50m

AND!  If possible, find the lower_bounds, upper_bounds, etc. without resorting to an external scripts or parsing gdalinfo outputs in the command line.


texts: osm2pgsql finishes processing without creating vertices_tmp table

texts:
I am using postgresql 9.1 database with posgis installed. I also have pgroute installed on same database. I am trying to use osm2pgsql tool to import data for routing. But osm2pgsql tool finishes processing without creating the vertices_tmp table, the output is like this:

Adding ways to database...
Creating topology...
NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
#########################
size of streets: 40234
size of splitted ways : 125002
finished

After i also try to do this step manually by executing this command:

SELECT assign_vertex_id('ways', 0.00001, 'the_geom', 'gid');

But it gives a similar error like this(with more details):

NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 15 at EXECUTE statement
ERROR:  query string argument of EXECUTE is null
CONTEXT:  PL/pgSQL function "assign_vertex_id" line 32 at EXECUTE statement

Any ideas how i can finish importing data for pgroute? I am using a city database(istanbul) osm file from http://metro.teczno.com .


texts: Is there an inbuilt arcpy method to insert a new line?

texts:
I need to extract a list of all the files and folders within a top level folder. My first thought was to use os.walk in a small python script as it is a simple task. However, I cannot run scripts on that folder as it is on the main corporate server.

I then took the same script, placed it in a Calculate Value inside a model and ran it successfully (the server allows all connections from ArcGIS). This still didn't give the result I wanted as I needed to catalog the spatial datasets as well.

I ended up using arcpy.da.Walk and wrote the list of files out to a text file. I'm having a problem with the new line character though, as was discussed in a fair amount of detail in this question. Every time I try to add it in, it gets stripped out immediately and an error pops up. Right now the code looks like this:

for filename in filenames:
  f.write(os.path.join(dirpath, filename))

It's writing each filename to the text file in "run-on" fashion, so it's quite useless. Is there some other way that I can insert a new line?

edit The issue at hand is forcing the inclusion of a newline character inside the model. I just happened to come across it while trying to create an "inventory". My question is not about how to get the list of files, but how to get the spacing to work correctly.


texts: Atmospheric correction of Landsat 7 slc-off images

texts:
I am interested in finding the change in land cover particularly forests.
I am working with slc-off images (no other sensor images available for my study). I have done the gap-filling in Erdas Imagine using focal-analysis and it worked fine. I have to further atmospherically correct the images using FLAASH in Envi.

What is the right sequence of process work?

Atmospheric correction before gap-filling or Gap-filling before atmospheric correction?


texts: What is the most efficient way to get the values of a previous row in a table using ArcPy?

texts:
I am very new to ArcGis and I would like to get a specific value of the previous row in a table.

I found this question here on gis.stackexchange.com and experimented a little bit. However, I found out this modified version:

Pre-Logic Script Code:

import arcpy

def getLastASTDate(tableName):
    return [r[0] for r in arcpy.da.SearchCursor(tableName, "date")][-1]

test =

getLastASTDate("name_of_table")

Unfortunately this version doesn’t work and is -- in my opinion -- very inefficient. Do you know any better solutions?


texts: Export shapefiles in QGIS to import to ArcMap

texts:
I have one question regarding export shapefiles from QGIS to import after to ArcMap 10.6. I export the shapefiles as Save data/ESRI shapefile. After when I try to upload in ArcMap this shapefile, I always got the same error: Could not upload the specific file.

The only difference between the shapefiles exported from QGIS is that instead of having .sbn and .sbx files, I have .qpj files.

Someone knows how to successfully exported in QGIS shapefiles for after importing to ArcMap?


texts: Link table values to ExtractByAttributes where_clause input

texts:
Using ArcGIS Desktop 10.7, I have 83 rasters that I need to loop through and perform an Extract By Attributes on using Python. Each raster has a range of values that the extracted pixels should be between. What I need to be able to do is insert the upper and lower bounding values from a table into the "where_clause".
ExtractByAttributes(in_raster, where_clause)
I'm not exactly sure how to link the appropriate values to the corresponding raster. Each raster is named after a "Map_ID" and the same "Map_ID" is located in the table. The table only has 3 columns, the Map_ID, upper and lower bounding values. The current python script I have is as follows (the 2 phrases "(Map_ID lower value)" and "(Map_ID upper value)" is what needs to call the values from a table):
# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

# Set environment settings
arcpy.env.workspace = "C:\data\Rast_conversion_Step_1"

# set output folder
outWorkspace = "C:\data\Rast_conversion_Step_2"

# create list of rasters to perform raster math on
RasterList = arcpy.ListRasters()

for item in RasterList:
      inSQLClause = "Value >= (Map_ID lower value) and Value <= (Map_ID upper value)"
      extractRast = ExtractByAttributes(item, inSQLClause)
      outputRaster = os.path.join(outWorkspace, item)
      ExtractByAttributes.save(outputRaster)


texts: Saving multiple raster sums to one Excel file on QGIS

texts:
I have a couple hundred raster files that I want to find the sum for in order to perform some further calculation. I want to have the sum of each raster file placed in one Excel sheet.

I used to work on ArcGIS, so usually I would just use Model Builder to iterate the procedure, then save it all to an attribute table and export to CSV.
But now I no longer have access to ArcGIS, and I realized that you can't iterate in the QGIS graphical modeler.
Is there a simple way to do this? I would prefer doing it without any scripting or code, because I just need to get this done quickly, but if it's absolutely necessary I don't mind.


texts: How to resolve AttributeError: 'list' object has no attribute 'listLayers' in Python 3 for ArcGIS Pro?

texts:
I have written several scripts for ArcGIS Pro and this error shows up occasionally and stumps me:

Traceback (most recent call last):
  File "<string>", line 10, in <module>
AttributeError: 'list' object has no attribute 'listLayers'

Simple code example:

aprx = arcpy.mp.ArcGISProject("CURRENT")
m = aprx.listMaps("pdf_template")
for lyr in m.listLayers():
    print(lyr.name)

line 10 is for lyr in m.listLayers():


texts: Python: how to speed up spatial search (nearest point)?

texts:
I'm working on a program that determines the closest location from a given point. The point cloud I'm testing against is very big (~800.000). However, this doesn't really explain why my implementation is so slow. This is my approach:

First, I created a spatial index for the point shape

pntshp.ExecuteSQL('CREATE SPATIAL INDEX ON %s' % table_name)

I defined an array of buffer distances to narrow down the search radius. Which of course also means that I have to create a buffer for each point (which might be expensive).

BUFFER_DISTANCES = ( 0.001, 0.005, 0.01, 0.02, 0.05 ) # 100m, 500m, 1km, 2km, 5km

Then, the buffer is used as a spatial filter

node_lyr.SetSpatialFilter(buff)

If the filter returns None the buffer distance will be increased.

for buffer_d in BUFFER_DISTANCES:
    buffr = get_buffer(xy_street,buffer_d)
    ...

Then I am calculating the distance to the points returned by the spatial filter

p=ogr.Geometry(ogr.wkbPoint)
p.AddPoint(xy[0],xy[1])

for feat in node_lyr:
    geom = feat.GetGeometryRef()
    d = p.Distance(geom)
    dist.append(d)

To get the closest point:

def get_closest_pnt(dist, node, how_close):
    mrg = zip(dist,node)
    mrg.sort(key=lambda t: t[0])
    try:
        return mrg[how_close]
    except IndexError, ierr:
        print '%s \ndist/node tuple contain %s' % (ierr,mrg)

It all works fine but is really slow. Creating a spatial index didn't show any effect, really. To calculate 100 points this implementations takes ~6,7 seconds. The program needs to be able to calculate the closest location for more than 2000 points as fast as possible. Any ideas on how to improve my approach?

EDIT

I tried different approaches to see where it gets me. I came across something very astonishing I want to share here.

I implemented a simple lookup algorithm as described here, and one of the solutions that where suggested (the sorted set approach).

The surprising fact is that performance is not only dependent on the implementation but even more so of the OSX. My original ogr/buffer algorithm turns out to be blazing fast on my OSX whereas it is painstaking slow on Linux (hence the question here).

Here are my results (100 runs).

Method       |     OSX        |  Linux Ubuntu
ogr buffer   | 0:00:01.434389 | 0:01:08.384309

sub string   | 0:00:19.714432 | 0:00:10.048649

sorted set   | 0:00:01.239999 | 0:00:00.600773

Specs Mac OSX
Processor 4x2.5 GHz
Memory 8 GB 1600 MHz

Specs Dell Linux Ubuntu
Processor 8x3.4GHz
Memory 7.8 GB

If someone can explain why these differences occur, please don't hesitate.


texts:
I am developing a plugin in QGIS using Python language. This plugin uses algorithms from Processing Toolbox and I save some results in the temporary folder assigning QgsProcessing.TEMPORARY_OUTPUT in the output field of each algorithm. However, my plugin is just working well when I return the temporary folder directory using QMessageBox. If I do not use this, my plugin does not work, so nothing is saved in the temporary folder. Why this happens?
Example of code working:
 clip_ndvi = Processing.runAlgorithm("saga:cliprasterwithpolygon", {
                                'INPUT': str(self.dlg.lineEdit.text()),
                                'POLYGONS': vl,
                                'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT})
                            out_ndvi = clip_ndvi['OUTPUT']
                            QMessageBox.about(self.dlg,'teste', str(out_ndvi))

Example of code not working:
clip_ndvi = Processing.runAlgorithm("saga:cliprasterwithpolygon", {
                            'INPUT': str(self.dlg.lineEdit.text()),
                            'POLYGONS': vl,
                            'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT})
                        out_ndvi = clip_ndvi['OUTPUT']

The difference is in the QMessageBox.