texts: OpenLayers3, display attribute information

texts:
I have OpenLayers OSM map with two WMS layers hotpoint and threat provided from GeoServer. I need to display attribute information from wms layers with popup or other appearing info table. I find this ask, but
don't understand this work. My code:

$(document).ready(function() {
view = new ol.View({
    center: [4701182.98765148, 7492051.764399836],
    zoom: 5,
    maxZoom: 18,
    minZoom: 2
});
var format = 'image/png';
var osm = new ol.layer.Tile({
    source: new ol.source.OSM(),
    visible: true,
    name: 'osm'
});
var mousePosition = new ol.control.MousePosition({
    coordinateFormat: ol.coordinate.createStringXY(2),
    projection: 'EPSG:4326',
    target: document.getElementById('myposition'),
    undefinedHTML: '&nbsp;'
});
var untiledhotpoint = new ol.layer.Image({
    source: new ol.source.ImageWMS({
      ratio: 1,
      url: 'http://192.168.255.197:8080/geoserver/geoportal/wms',
      params: {'FORMAT': format,
              'VERSION': '1.1.1',
               LAYERS: 'geoportal:hotpoint',
               STYLES: '',
      }
    })
});
var hotpoint = new ol.layer.Tile({
    visible: false,
    source: new ol.source.TileWMS({
      url: 'http://192.168.255.197:8080/geoserver/geoportal/wms',
      params: {'FORMAT': format,
              'VERSION': '1.1.1',
               tiled: true,
               LAYERS: 'geoportal:hotpoint',
               STYLES: '',
      }
    })
});
  var untiledthreat = new ol.layer.Image({
    source: new ol.source.ImageWMS({
      ratio: 1,
      url: 'http://192.168.255.197:8080/geoserver/geoportal/wms',
      params: {'FORMAT': format,
              'VERSION': '1.1.1',
              LAYERS: 'geoportal:threat',
              STYLES: '',
      }
    })
 });
  var threat = new ol.layer.Tile({
    visible: false,
    source: new ol.source.TileWMS({
      url: 'http://192.168.255.197:8080/geoserver/geoportal/wms',
      params: {'FORMAT': format,
               'VERSION': '1.1.1',
               tiled: true,
            LAYERS: 'geoportal:threat',
            STYLES: '',
      }
    })
  });
var map = new ol.Map({
    target: 'map',
    controls: ol.control.defaults().extend([
        new ol.control.ScaleLine(),
        new ol.control.ZoomSlider()
        ]),
    layers: [osm, untiledhotpoint, hotpoint, threat, untiledthreat],
    view: view,
    renderer: 'canvas'
});
    map.addControl(mousePosition);
    var attributeData = function(pixel) {
        var feature = map.forEachFeatureAtPixel(pixel,  function(feature, layer){
            return feature;
        }, null, function(layer) {
            return layer === hotpoint;
        });
        var info = document.getElementById('info');
        if (feature){
        info.innerHTML = '<div>' + feature.get('point_id')+'</div>'
    } else {
        info.innerHTML = '&nbsp;';
    }
};
});


texts: Load single image by map API and draw overlays

texts:
I have an image which  I want to make it zoomable and dragable in browser, and draw some overlays like lines and markers on it.

So I wonder if I can use the map api like Google map or openlayers to load the image. Then I can use the draw feature provided by the API to add overlays.

Since the image size is not fixed ,some of them are almost 100M. So I may have to generate tiles first, but the image have no coordinates, so I do not know if it is possible.

Any suggestions?


texts: Finding all line segments lying inside polygon and what % of line inside polygon using Google Fusion Tables?

texts:
I was practicing with fusion tables sample code and while playing with some dataset, I was thinking how to achieve this. So if you see attached image.

Red blocks : These are SLA code or SLA areas. And if you click on it, will tell you total patients in that area having (COPD , Heart Disease, Lung cancer).

Yellow  lines : Boundary between different SLA regions.

Blue lines: Vic roads. When you click any line inside the red box. It will tell you road name, total vehicles.

Now, I want all the blue lines inside each red blocks basically I want total vehicles which is sum of all the vehicles on each blue line. I also want to find what % of whole blue line is lying in each red .


texts: GP service (from Python script tool) time-out: how to execute finally: statement

texts:
I am trying to build more fail-management logic into some of the GP services I develop. As a part of the business logic within the script tool I build, I write some data into the tables and thus lock access to certain parts of the geodatabase to let the currently running GP task to be executed undisturbed.

However, in case there are some issues which will force the GP service instance or ArcGIS Server Windows service to crash (something that is out of my hands), I want to be able to perform some clean-up work (release locks, clean-up scratch geodatabases, write to tables, etc). This is where I usually use finally: statement when developing Python script tools for use in ArcGIS Desktop.

However, the GP service does not seem to reach the finally: if it is timed-out or crashed.
The Python script tool source code is below:

import arcpy
import time
try:
    arcpy.AddMessage("Hello 1 second")
    time.sleep(10)
    arcpy.AddMessage("Hello after 10 sec")
except:
    pass
finally:
    arcpy.AddMessage("GP service timed out!")

I have run this tool in ArcMap:

Trying to run the GP service task (published from the tool result) that of GP service with time-out time for 9 seconds:

So the GP service is not able to reach finally statement and print "Hello after 10 sec". If this is not possible to implement within the finally: statement in Python script tool source code, how would you approach this otherwise?


texts: PostGIS- GROUP BY for continuous data values?

texts:
In the sample code for the PostGIS ST_ConvexHull function, multiple polygons would be created based upon different, discrete values found in the field "disease_type".

--Get estimate of infected area based on point observations
SELECT d.disease_type,
    ST_ConvexHull(ST_Collect(d.the_geom)) As the_geom
    FROM disease_obs As d
    GROUP BY d.disease_type;

Lets say there was a numeric field (lets call it "myfield") with continuous values from 0-5.  How would you produce similar results as GROUP BY but defining your own breakpoints (for example: 0-1.25, 1.25-3.5, 3.5-5)?


texts: Using ESRI geoprocessing services from a non-ESRI client (Python)?

texts:
I'm thinking it should be possible to consume an ESRI Geoprocessing service from a python script without invoking ArcPy or needing an ESRI installation on the client machine.  For example I create a geoprocessing service that takes a series of point location and returns a csv file with the values under them as a file.  I would like to include this service call in a python client application that will be running on machines that don't have any ESRI installation.  Couldn't you create the rest call to kick off the service and consume the result without arcpy?

Does anyone have an example of how you might do this?

Thanks,
Colin


texts: How to integrate eclipse (PyDev) with Python Add-In development (ArcGIS 10.1)?

texts:
In ArcGIS 10.1, Python Add-Ins make various buttons, toolbars, extensions, etc, available to those already familiar with Python and OO in general. I am in the process of developing a Python Add-In and normally I would use eclipse with PyDev to develop scripts or small python applications. However, with the prescribed/predefined Python Add-In workflow (Python Add-In Wizard), I am finding it difficult integrating my eclipse workflow.

The Python Add-In relies on a script that packages the add-in according to a predefined config.xml and folder/file structure, seen below:

Python Add-In\
    Images\
        yourimages.png
    Install\
        PythonAddIn_addin.py
    config.xml
    makeaddin.py

How would you integrate eclipse PyDev without interfering with the python add-in wizard and prescribed development folder structure?


texts: GDAL - Specify span of elevation values for conversion to Terragen format?

texts:
I'm attempting to use the gdal tools to convert some geotiff files I have to the Terragen format. I was hoping to use either gdal_translate or gdalwarp in a similar fashion to this post. (though I should admit I am not sure which tool may be more appropriate for my task)

However, when I issue either command it exits with an error, in fact, the same error:

$ gdal_translate -of Terragen src.tif dest.ter
Input file size is 8743, 6159
0ERROR 1: Inverted, flat, or unspecified span for Terragen file.

$ gdalwarp -of Terragen src.tif dest.ter
Creating output file that is 8743P x 6159L.
ERROR 1: Inverted, flat, or unspecified span for Terragen file.

I quickly stumbled upon a post (which also points to gdal documentation) stating that I need to specify the span of elevation values:

  Use the Create call. Set the MINUSERPIXELVALUE option (a float) to the
  lowest elevation of your elevation data, and MAXUSERPIXELVALUE to the
  highest. The units must match the elevation units you will give to
  band::SetUnitType().

Were I developing my own utility that used libgdal, this information would be applicable.  However, instead of developing a new utility to achieve the conversion - I'd rather use the existing gdal tools if possible ( gdal-translate, gdalwarp).

Is it possible to specify the min/max elevation values using any existing gdal tools?


texts: How do you restrict road access by time in ArcGIS Network Analyst?

texts:
I was wondering if there's a way to restrict road access based on time of day (i.e., create a time restricted road).
I'm familiar with restricting roads based on vehicle weight/height, but time restriction is a little different-- there's no time unit in ArcGIS (or at least I'm not aware of it).

I was thinking of using python to script the restrictions. I'd create 2 new fields (double) in the road network attribute table, 1 for when the time restriction starts and 1 for when it ends (restrict_tstart, restrict_tend).
For example if the road is restricted between 10am - 2:30pm, I would put 10 in restrict_tstart and 14.5 in restrict_tend.
Then I'd write a python code that returns True when the real time falls between restrict_start and restrict_tend.

Does that work? I'm not even sure where to put that script in ArcMap.

Suggestions, anyone?


texts: Using rainfall data (raster) as an input to determine stream lines in Arc Hydro?

texts:
I was wondering if it is possible to use rainfall data (raster) as an input to Arc Hydro for determining stream lines?

I am learning Arc Hydo and have gone through a couple of tutorials. As far as I can see, the tutorials mostly talk about DEM editing (selecting/filling sinks, AgreeDEM etc.). But I can't find a place where rainfall data can be used to determine the amount of water entering leaving cells. It seems as if the stream definition were only based on topography and the number of upstream cells.

My questions are:

Is there a way (and documentation/tutorial) to use rainfall and/or landcover (e.g. soil vs concrete) data as an input to Arc Hydro for determining streams and watersheds?

If yes, how do modelers use such input? (e.g. by defining some weight values for each cell.)

If not, are there other ArcGIS or open source package similar to ArcHydro but can take into account precipitation?

P.S. I searched on GIS SE, and there is related question about Calculating watershed runoff with landuse and rainfall data? how to compute rainfall. My question concerns how to use rainfall data.


texts: Rasterise polygon layer using QGIS 3.0.3

texts:
I am trying to rasterise a simple polygon layer this file using QGIS but cannot get a sensible output.

The documentation doesn't appear to have been updated for version 3 yet, and does not match.

The options I have are
.

What options I be using to rasterise the polygon layer by attribute LandID?

(I got the result I expected swapping to 2.8 - but I'd like to figure out how to do this with the new version)


texts:
I've to compare two polygon layers, 'land_use' and 'land_parcels', in PostGIS. My aim is to create linestrings, where boundaries of 'land_use' polygons do not match boundaries of 'land_parcel' polygons.

A similar question has already been asked as How to create lines to visualize differences between polygon features in PostGIS?.
Can anyone help me adapting the query posted by Mike T?
SELECT row_number() OVER () rn, *
FROM (
  SELECT (ST_Dump(ST_SymDifference(ST_Snap(luse, lparcels, tol), ST_Snap(lparcels, luse, tol)))).*
  FROM (
    SELECT ST_Union(DISTINCT land_use_only.geom) luse, ST_Union(DISTINCT land_parcels_only.geom) lparcels, 1e-5 tol
    FROM (
      SELECT ST_Boundary(geom) geom FROM land_use
      EXCEPT SELECT ST_Boundary(geom) geom FROM land_parcels
    ) land_use_only,
    (
      SELECT ST_Boundary(geom) geom FROM land_parcels
      EXCEPT SELECT ST_Boundary(geom) geom FROM land_use
    ) land_parcels_only
  ) s
) s;


texts: Two maps, one with classes and one with no classes. get the sum of one map into the classes of the other

texts:
I have two maps. One map has some classes and the other one is a raster with no classes.
To illustrate, the map with the classes: 0-1, 1-2 and 3-4! The other map is a raster with values from 0-1000.

The maps are overlaying. I want to find out the sum of the second map within the classes of the first map, i.e. what is the sum in the first class 0-1 of the values of the second map?

ArcGIS 10.2


texts:
I'm  new to web mapping, but I've managed to put together a fairly simple web map using Leaflet:
http://www.hcswgis.com/globalmapper/

It has most of the functionality that I want using the available plugins; One feature that it's lacking is a way to search by a parcel ID, which in turn puts a marker where that parcel is, and zooms into it.

In my web map there is a parcel Mapbox layer, but I don't believe there is a way to do what I'm trying to accomplish with that layer because it's tiles and a UTFGrid.

I thought about adding the parcel layer onto my GeoServer. From that I thought I could potentially use the plugin https://github.com/stefanocudini/leaflet-search to query the GeoServer parcel layer to find whatever parcel ID is being searched and zoom into it.

Is that possible?

I'm open to any alternative solutions that accomplishes the goal. Keep in mind we're talking about around 80,000 features to search.


texts: Why does Google Earth Pro change my raster resolution?

texts:
I have noticed this problem occurring when I try to import rasters into Google Earth Pro from ArcGIS Desktop 10.4.1. I can take a 1/3 arc second resolution DEM and it looks great in ArcGIS until I use the layer to KML in which case it changes the resolution to something strange.

When I measure the pixels in Google Earth, it's about 150 meters per pixel side where it's only 10 meters in ArcGIS. I have played around with exporting at a higher DPI and a 2048 size of returning image. Anyone have any ideas on how to get a DEM into Google Earth in its native resolution?


texts:
I have three geo-referenced images of a museum layout that I have been able to offset by "x" number of meters in ArcScene and ArcGIS Pro. I want to publish the three images online and have it look like below:
I am not able to find a way to add offset to the layers in ArcGIS Online Scenes. I have tried to exported directly to a Exporting ArcScene Documents to 3D Web Scenes and viewed in Esri CityEngine Web Viewer. However the image quality is very low and transparent image backgrounds are not allowed.
Unfortunately this does not allow for easy navigation and will not work for us. I want to do this in ArcGIS Online if possible but I am open to other solutions that would allow other GIS features to be added to the images.


texts: Deleting polygons within polygons using Model Builder

texts:
I have smaller polygons -within bigger polygons-, and they supposed to by "holes" in the bigger ones. Based on the answers of this question I know how to identify them, and also how to get rid of the small polygons filling up the holes.
I also managed to make a model based on this. My problem is that it only works if the source feature classes are added to the table of contents and I put them in the model manually. Is there an other way to do this process as the "inside" part of a way bigger model?
Here is my model for easier visualization of my problem:

The source and the source copy are identical
Select Layer By Location: I select the features that are comletely within the other features, which results the "inner rings" of the polygons.
Select Layer By Location(2): I simply switch the selection, so I can remove the unnecessary features, which happens in Delete Features.
Erase: I use the previous results as the erase feature with the copy of the original feature class as the input feature.

Now this works perfectly as I said, but only if the source and source copy are already added to the Table Of Contents, and only if I set them in the model manually every time.
Is there a way to do the same process within a model, with feature classes as the results of some other tools?
I am using ArcGIS Desktop 10.8.


texts:
I have a feature created from Well Known Text as following:

var feature = format.readFeature(wkt , {
    dataProjection: 'EPSG:25830',
    featureProjection: 'EPSG:25830'
});

When I add it to a layer and try to center in the map like this:

this.vectorSource.addFeatures(feature);
this.map.getView().fit(feature.values_.geometry.extent_, {duration: 2500});

I get an AssertionError. When analyzing the feature I find the extent is set to infinity:

extent_: (4) [Infinity, Infinity, -Infinity, -Infinity]

  How can I calculate the feature extent? How can I zoom and center the
  created feature on the map?


texts: Two different ways of labeling in QGIS Layouter

texts:
I've made a map with points and labels related to those points. In one area, the points are very dense, therefore, I'd prefer to create one small scale map and in the corner a large scale map, to get more detail.

Now, my problem is that - for the dense part - labels should only be put on the large scale map, for the other points, the label should be present on the small scale map.

How would you solve this?

I can choose labels for one map (like in: How to only show labels for an arbitrary selection of items?), make a layouter for the map, and subsequently do the same for the other map. But than, if you want to change something to the first, labels are updated again to the actual state of your qgis project.


texts:
I am trying to change the DPI of QGIS3 (3.6) in the layout manager.
However, as the Figure below shows (lower left corner), the DPI (300 by default) seems to be read only.

I tried to click the textbox. But it's not allowing input.
There is no spinner control QGIS 2 either as said about the the older QGIS2 in GIS.SE:

Adjusting QGIS Print Composer output resolution?

How can change the DPI from the QGIS 3 GUI?


texts: ArcGIS Extract by Attributes Error 010411

texts:
I am attempting to process a raster into a polygon, with the extent of the polygon being where the raster contains prescribed allowable values. I have done this successfully before using the following steps:

Creating a raster of 1 for allowed values (e.g. VALUE >= 1 & VALUE <= 5) and 0 for not allowed values using the Raster Calculator tool.
Use the Extract by Attributes tool with the setting Value = 1 creating a raster of allowed values only.
Converting this new smaller raster to a polygon using the in-built ArcGIS Raster to Polygon tool.

I have hit an error I have not before in this process. I am starting with this raster, where I have use Classified symbology to show allowable and not allowable values:

I then run Raster Calculator with the following settings:
Con("myraster.tif">=685,1,NULL)  &  Con("myraster.tif"<=735,1,NULL)

Output raster: "t_one_WA"

This produces the expected result:

I then attempt to run Extract by Attributes with the following settings:
Input raster: t_one_WA
Where clause: VALUE = 1
Outputraster: mylocaldrive_t_two_WA

This returns a number of errors, starting with Error 010411 "Error in accessing symbol __1000000 in the symbol table." Looking at the "Help" guide for this error suggests exporting the raster to a local drive, which I am and am not having success. The full "Messages" of Extract by Attributes is here:

(1) What are the possible causes of this problem?
(2) Is there a different way to achieve my desired result from my initial raster?


texts: Add a category variable that gives a unique id for each segment of 10m of each street

texts:
I have a big csv (12Go) that contains this informations of trip data for each point in 1s interval: lon, lat, timestamp and street_id.
My aim is to segment street to 10m long segments and assign a unique segment id to each point based on its position.
Here is an example of what my data looks like using df.head():
    timestamp           Lat       Lon       street_id
0   4/1/2014 0:11:00    40.7690   -73.9549  140
1   4/1/2014 0:17:00    40.7267   -74.0345  50

I'm using goepandas as follows to get a spatial dataframe:
from geopandas import GeoDataFrame
from shapely.geometry import Point

geometry = [Point(xy) for xy in zip(df.Lon, df.Lat)]
df = df.drop(['Lon', 'Lat'], axis=1)
gdf = GeoDataFrame(df, crs="EPSG:4326", geometry=geometry)

I just dont know how to begin, I know I have to group by street_id then calculate the distance between each point and the begining of the street, but Ive failed to find how to find how to do so
The expected result should look like:
    timestamp           Lat       Lon       street_id segment_id
0   4/1/2014 0:11:00    40.7690   -73.9549  140          5
1   4/1/2014 0:17:00    40.7267   -74.0345  50           1

Segment_id should be a unique identifier for each street segment of 10m long


texts: How do you view data contained inside hidden directories with QGIS browser?

texts: I am attempting to browse data in a hidden directory, but the QGIS Browser tree view doesn't show in hidden folders? I know, I know...it's hidden. But is there a command-line switch to override this behavior? The folder in question is APPDATA, so I don't want to make it visible.


texts:
Does anybody know how to perform a full outer join with ArcMap?

I have the following two tables:

Table A
Joinkey
1
2
3

Table B
Joinkey Other_field
1 1
1 2
1 3

The result should be something like:

Joinkey_A Joinkey_B Other_field
1 1 1
1 1 2
1 1 3

But ArcMap does the following (inner join):

Joinkey_A Joinkey_B Other_field
1 1 1
2 2 1
3 3 1

In my opinion this means, I need a full outer join but I cannot find any solution for this.


texts: ST_contains between two tables returns all results of one table multiplied by the number of results from other table

texts:
I have one PostGIS table of polygons like this:

polygon_id | polygon_geom

And another of points like this:

point_id | polygon_id | point_geom

I want to find all the points from the points table that are within a specific polygon from the polygons table. I tried this query:

SELECT points.point_id
FROM   points,
       polygons
WHERE  St_contains((SELECT polygons.polygon_geom
                    FROM   polygons
                    WHERE  polygons.polygon_id = 2905), points.point_geom) = 't';

I expect 21 point_id's, but instead it gives me the entire polygon table with each line repeated 21 times. How can I get just the 21 rows from the points table? And furthermore, how can I then assign the polygon_id that these points are within to the points.polygon_id column for those 21 rows in the points table?

EDIT: If I use this query to check if one polygon is in another from the polygon table, this query works perfectly. So I feel like I'm close but fundamentally misunderstanding something.

EDIT 2: changed polygon.id column to polygon.polygon_id


texts:
I want to do some isochrones with pgrouting, and I am failing in setting up the network in the first place.

osm2pgrouting doesn't find my postGIS installation. Also, the error message is not very helpful (ERROR: postGIS not found) and I don't know what to do (reinstalling already tried).

What is your advice?

osm2pgrouting -f lebanon-latest.osm -d routing_lb -U jb -W
***************************************************
           COMMAND LINE CONFIGURATION             *
***************************************************
Filename = lebanon-latest.osm
Configuration file = /usr/share/osm2pgrouting/mapconfig.xml
host = localhost
port = 5432
dbname = routing_lb
username = jb
password =
schema=
prefix =
suffix =
Don't drop tables
Don't add nodes
***************************************************
Connecting to the database
host=localhost user=jb dbname=routing_lb port=5432 password=
connection success
ERROR: postGIS not found

In PostgreSQL, the postGIS extension is installed:

postgres=# create extension postgis;
CREATE EXTENSION
postgres=# create extension pgrouting;
CREATE EXTENSION

Test:

postgres=# create extension postGIS;
ERROR:  extension "postgis" already exists

Not surprisingly, no tables in the database have been created:

postgres=# \c routing_lb
You are now connected to database "routing_lb" as user "postgres".
routing_lb=# \dt
No relations found.


texts: ArcGIS JavaScript v4 - KML via FeatureLayers

texts:
Now i'm aware there's no KML module in 4.0 beta but this blog http://odoe.net/blog/diy-kml-in-arcgis-js-api-4-beta/ has a pretty nice way to overcome this. Unfortunately not working for me.

[fiddle] https://jsfiddle.net/ytsycg5r/

1) I use "//utility.arcgis.com/sharing/kml" to convert kml to feature set
2) use GroupLayers to create featurelayers and add to map.

I get error "esri/request: the first parameter should be a url string"


texts:
Sorry for my lack of knowledge.

I have been using GIS for a while but just for analysis and don't know much about the data side of the business.

Recently my boss asked me to think about setting up a GIS platform that should be able to handle this type of data appended to property location:

Co-ordinate data
Location value
Basic hazard data
Model output data
Detailed reports from our Engineering field force

Potentially the number of properties could runs into millions and it could be interesting if the engeneers on the field could have access to the maps and input information live into the database.

Also maps should be made available to colleagues and stakeholders through an online platform.

I have always been using ESRI but never used ArcGIS Server and it would be great if somebody with more experience than me could put me on the right path.

How would you set up the system?

Are there openesource alternatives?


texts:
So, I'm writing an web mapping application using OpenLayers in WebForms in .NET 4.5.

It was working wonderfully until my browser upgraded to IE9. Now the JavaScript fails. After researching I found that the problem seems to be that IE9 changed how it handles certain kinds of JavaScript. So the solution, until the community can rewrite the JavaScript libraries, seems to be to have the browser run in compatibility mode. Unfortunately, I can't seem to get compatibility mode to work correctly.

Currently my application has one page and this is the header:

<%@ Page Language="C#" AutoEventWireup="true" CodeBehind="default.aspx.cs" Inherits="AgPlaceHolder._default" %>

<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">
<head runat="server">
<meta http-equiv="X-UA-Compatible" content="IE=EmulateIE8" />
<title>Agrilogic Place Holder</title>
<script src="OpenLayers.js"></script>
<script src="AgInit.js"></script>
<script src="AgLabels.js"></script>

<style>
    html, body
    {
        width: 100%;
        height: 100%;
        margin: 0;
        padding: 0;
    }
</style>
</head>

It didn't work. So I found this alternative to add to the web.config file:

<system.webServer>
  <httpProtocol>
    <customHeaders>
      <add name="X-UA-Compatible" value="IE=EmulateIE8"/>
    </customHeaders>
  </httpProtocol>
</system.webServer>

This also didn't work.

I've tried some variations I've seen used around, but none of them seem to be working. Does anyone see what I might have done wrong? Or has another solution to this problem?

Thanks.

On Edit:

I'm using OpenLayers 2.12.

The error is:

Unhandled exception at line 670, column 407 in http://localhost:57639/OpenLayers.js

0x800a138f - Microsoft JScript runtime error: Unable to get value of the property 'resources': object is null or undefined

The code for setting up the map is:

map = new OpenLayers.Map({
    div: "map",
    layers: layerArray,
    controls: [
        new OpenLayers.Control.LayerSwitcher({ roundedCorner: true }),
        new OpenLayers.Control.PanZoomBar(),
        new OpenLayers.Control.Navigation(),
        new OpenLayers.Control.ScaleLine({ geodesic: true }),
        new OpenLayers.Control.KeyboardDefaults(),
        new OpenLayers.Control.MousePosition({
            displayProjection: new OpenLayers.Projection("EPSG:4326"),
            numDigits: 5,
            emptyString: '',
            prefix: "x: ",
            separator: " y: "
        }),
        //new OpenLayers.Control.EditingToolbar(editableLayer),
        //new OpenLayers.Control.SelectFeature(editableLayer, { toggle: true, clickout: true }),
        //new OpenLayers.Control.DragFeature(editableLayer),
        overviewMapControl
    ]
});

The java script works fine, except for another issue I have posted a question on over in GISStackExchange, but that's not preventing the maps from showing. The problem there is I can't drag the label-points the program creates. But now with the upgrade of the browser I can't even display the maps anymore. The page just fails.

The newly built page has the same errors. When I remove Bing:

0x800a138f - Microsoft JScript runtime error: Unable to get value of the property 'wrapDateLine': object is null or undefined

When I include Bing:

0x800a138f - Microsoft JScript runtime error: Unable to get value of the property 'resources': object is null or undefined

I am at the mercy of your expertise.


texts: Merge multiple layers in different folder

texts:
I am working with Korean topo map layers.
They are 1:1000 scale and divided by each map section.

they are also stored in different folder.

Since a layer shows road has same name (e.g., a0010000_polygon.shp), i added the district name as a prefix (e.g., 376120547_blah blah_..._a0010000_polygon.shp)

Here is the problem.
now, i need to merge each layers that have same suffix (i.e., merging all layers that ends with a0010000_polygon.shp in each folder) to make one large big huge layer

What i did is, iterate thru all the folder and using wildcard and make feature layer (add to the TOC). Then, drag the layers to merge toolbox..

there are at least 180 layers in each folder, and doing 180 times of this is kinda crazy.

Are there better way to solve this issue?

I came up with this model ... (using collect value) but still i need to adjust wildcard to iterate thru layers have same suffix :(

Oct. 16, 2016 I add this comment

my final (?) model looks like this. (it ain't final until it really works)
what i did is that i create a table that contains all possible layers (See 'library' field)

then i iterated thru the table row and get the row value and put that in my model.. so that it merges all layers that have same name

the problem is that the shape layers in a folder doesn't have all layers in the library table (e.g., the library has all possible layers, which is 108. but shape layers in a folder has only 40 or sometimes 60)
and therefore, it returns merging error lol

so i was thinking about using try / except python script and put that in a model
check this --> [http://resources.esri.com/help/9.3/arcgisdesktop/com/gp_toolref/automating_your_work_with_models/branching_colon_implementing_if_then_else_logic.htm][5]

spent almost a week to solve this issue with modelbuilder and it has some limit i think.
I will come up with better idea if i can. (so far the best way / easiest way is insert arcpy script try/except)
but if there is anybody who knows more precise and elegance way to solve this, please share with us.

after i tried 'dru87''s code, i finished the work.
I typed line by line into the console and run and it only returned the last string in the list.

instead, you must copy and paste the entire code and run it (of course, you need to change the workspace and stringSearch)

awesome. and nice. still i have the modelbuilder here and i prefer modelbuilder because i suck at python and if anybody can improve the model, please let me know!


texts: Georeference an image in WGS84 where I have the center lat/lon and image dimensions

texts:
I have a drone image which I have warped into the proper orientation and shape and now I'm trying to georeference this image. My image pixels are 0.2 meters, my image dimensions are 750x750 pixels, the center coordinate is 34.7150963, -76.6732313, and the x axis runs west to east and y axis runs north to south. I've provided code on my attempt but it may well be the wrong approach and really I just need to georeference this image based on the above parameters.

I've tried many approaches with rasterio such as trying to estimate the corner positions based on the offset from the center and then using the rasterio.transform.from_bounds() function. I found the bounds by estimating the offset in lat/lon based on 750 pixels*0.2m then divided by half to get the offset in meters from the center. I then tried to calculate lat/lon offsets from these offsets in meters but it seems to lead to a lot of error. There has got to be a simpler way to just georeference an image based on a center coordinate when it is already projected for WGS84.
First to get the coordinates of the corners of the image bases on a conversion from meters to lat/lon:
import math
scaling = 0.2 # pix size in meters

def calc_lat_lon(lat, lon, dn, de):
    #Earth’s radius, sphere
    R=6378137

    #Coordinate offsets in radians
    dLat = dn/R
    dLon = de/(R*math.cos(math.pi*lat/180))

    #OffsetPosition, decimal degrees
    lat1 = lat + dLat * 180/math.pi
    lon1 = lon + dLon * 180/math.pi
    return(lat1,lon1,0)

def output_corners(camera_lat, camera_lon, offset):
    #Position, decimal degrees
    lat = img.latitude
    lon = img.longitude

    top_left = calc_lat_lon(lat, lon, offset, -offset) # top left corner
    top_right = calc_lat_lon(lat, lon, offset, offset) # top right corner
    bottom_right = calc_lat_lon(lat, lon, -offset, offset) # bottom right corner
    bottom_left = calc_lat_lon(lat, lon, -offset, -offset) # bottom left corner

    return(np.array(top_left, top_right, bottom_right, bottom_left))

coords = output_corners(img.latitude, img.longitude, img.shape[0] * scaling / 2)

coords returns:
array([[ 34.71577004, -76.67405094,   0.        ],
       [ 34.71577004, -76.67241166,   0.        ],
       [ 34.71442256, -76.67241166,   0.        ],
       [ 34.71442256, -76.67405094,   0.        ]])

Then to do the actual transform and try to georeference:
n = np.max(coords[:,0]) # max latitude aka north
s = np.min(coords[:,0]) # min lat aka south
w = np.min(coords[:,1]) # max lon aka west
e = np.max(coords[:,1]) # min lon aka east

import rasterio
from rasterio.warp import reproject
from rasterio.enums import Resampling

with rasterio.Env():

    # open the original image to get some of the basic metadata
    with rasterio.open(path_name, 'r') as src:
        profile = src.profile

        dst_crs = "EPSG:4326"

        tsfm = rasterio.transform.from_bounds(w, s, e, n, 750, 750)

        profile.update(
            dtype=rasterio.uint16,
            transform = tsfm,
            crs=dst_crs,
            width=img.shape[0],
            height=img.shape[1]
        )

        print(profile)

        with rasterio.open('new_georef_img.tif', 'w', **profile) as dst:
            dst.write(img.astype(rasterio.uint16), 1)

But again this leads to a fair amount of error and when I do this transform with GCPs or with other methods they all differ and lead to error so I'm looking for help in a better approach.


texts:
I am trying to create an autofill format on QGIS, based on an atlas.

There will be some dynamic text fields and some boxes which should be checked with an "x" in case the atlas geometry intersects a group of other geometries. Something like "IF $atlasgeometry intersects layer X OR Y OR.. OR Z, THEN check the box". Those layers are grouped, so it would be nice if i could call them all together in the script. The point is that i cannot find a single tutorial explaining how i should call map layers within the formulas, all the examples i see are about single geometries, and i have hundreds of them, so that would not work.

Could anyone explain me how that should work? And how the formula stucture should be?

Let me try to explain:

We have a map layer on which there are some objects stored: they describe the geographic positions of some buildings. Each geometry is a building. Let's call the layer "Projects"
We have some other layers containing the perimeters (=polygons) of the areas where some kind of jobs cannot be done or must have some features. Let's call them "Restraints A" "Restraints B" ...
I need a script which says "IF "projects" intersects Restraints A, then write X in a certain text box.

The goal is creating a format which describes in a single sheet all the restraints i must respect in the definition of the project.


texts:
I am trying to calculate the distance using the ST_Distance function in the postgres database.

I am using following query:

SELECT ST_Distance(
  ST_GeogFromText('SRID=4326;POINT(52.3551151038034 4.87660185549649)')
  , ST_GeogFromText('SRID=4326;POINT(52.3549985339677 4.87622001617791)')
);

For this the distance is returned as 44.159975697 meters.

When I am comparing this value with other sources, for example Haversine formula explained in this post.
The Haversine formula return the value as 0.0289913126501404, which is roughly 29 meters.

What am I doing wrong with my query here?


texts:
Where do I find this plugin?  ( I am on v. 1.6 - it does not seem to be installed and I could not find it with the python plugin installer).
-Robert


texts: How to manually change from GPS (WGS84) to web based maps projection systems?

texts:
I have collected location point data and wish to display this on a web based map. It was collected using a Garmin Etrex GPS using the WGS84 datum. I wish to display the points using Ushahidi and have done so successfully. The site is here.

http://www.flamingocensus.com/ushahidi/

Unfortunately the points do not correlate to what I know they should be, I believe it is a projection error. The points are supposed to be on the edge of the lake. Any suggestions or pointers on how to translate the points in the current projection to my new one.


texts:
In QGIS 2.18.16

I have 2 layers, one is for points with information about an accident and the other layer is a line string that represent "the road that I take care off".

I want only to extract or filter the information of points of an accident (yellow dots) that happen in the road that I take care off (blue lines).

I try with the "Joint attributes by location" tool but it doesn't work.

Please help!


texts: Get eigenvalues of large point cloud using lidR

texts:
I have some very big TLS point clouds (up to 400 million points) and want to assign the eigenvalues to each point in order to later calculate geometric features like planarity or sphericity. I use R and the package lidR following lidR book guidances. I realized that the function point_metrics() in combination with fast_eigen_values() is too slow. It takes ~15 minutes to compute
# load data
cloud_raw <- readLAS(path_points)

# because eigen() is really slow, use C++
Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::vec eigen_values(arma::mat A) {
arma::mat coeff, score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
return(latent);
}")

metric_geometry_features <- function(x, y, z) {
  xyz <- cbind(x, y, z)
  cov_matrix <- cov(xyz)
  eigen_matrix <- eigen_values(cov_matrix)
  geometries = list(
    planarity = (eigen_matrix[2] - eigen_matrix[3]) / eigen_matrix[1],
    linearity = (eigen_matrix[1] - eigen_matrix[2]) / eigen_matrix[1]
  )
  return(geometries)
}

metrics <- point_metrics(cloud_raw, ~metric_geometry_features(X,Y,Z), k = 20)

The function segment_shape() is way faster and also computes eigenvalues in the background. However, looking at the C++ Code behind lidR - although I can't write C++ - it seems to be that it is written in a way that it only able to return booleans, not doubles.
Does someone know how to solve this problem? Was there already a version somewhere implemented, which enables me to use the spatial indexing to return eigenvalues without needing to learn C++ just for this? My next resort would be to try using voxel_metrics(). However, I would prefer to have one set of eigenvalues for each point.


texts: "User memory limit exceeded" when selecting a single feature in a feature collection?

texts:
I have a feature collection of 18000+ polygons. I would like to get the pixel values of each pixel in every polygon that intersects Landsat imagery. The problem I am currently having is when I try to access one feature in the feature collection, I am getting the "user memory limited exceeded" error. Is this feature collection too much for Earth Engine to handle?

var table = ee.FeatureCollection("users/adrianom/SN_Meadows");
var size = table.size();
var fc_list = table.toList(size);
var feature = fc_list.get(0);
print(feature);


texts: QGIS Modeler: Additional Option don't work

texts:
I'm trying to automize clipping a raster obtained from an online source (WMS). It's all ok, when I do it step by step, but I need the final raster file to be compressed; and that's the point that doesnt work in the Modeler.
That's what I do:

creating a buffer (1500m) around my AOI (vector polygon)
"downloading" the WMS raster by using "Convert map to raster" (extent: buffer #1)
creating a buffer (1000m) around my AOI (vector polygon)
clipping raster (#2) to the second buffer (#3) by using "Clip raster by mask layer": compared to "Convert map to raster" I can
choose between different compression settings (section "additional
creation options") -> "low compression"

These 4 steps I put in the Modeler. It all works like a charm, but the final raster is not compressed :-( In the dialog for "Clip raster by mask layer" I can't choose any compression presets. Instead, two fields "additional options" and "additional command-line parameters" are avaible, when "Show advanced parameter" is activated at the top of the dialog. I figured out that it should be somthing with "COMPRESS=PACKBITS".
Here's what I tried in both lines without an effect on the size of the output raster:
-co COMPRESS=PACKBITS

COMPRESS=PACKBITS

'COMPRESS=PACKBITS'

gdalwarp -co COMPRESS=PACKBITS

gdal_translate -co COMPRESS=PACKBITS

It seems, that all these options are ignored. So I assume, the ways I tried to use them are wrong.
I'm looking foward to hints and explanations.


texts: When saving a .qgs file in QGIS, another appears. Is this a backup?

texts:
When saving a project file (.qgs), which when executed is automatically loaded with qgis.bat, there's another .qgs file which appears in the same directory of the original project file, although this one is not loaded with the qgis.bat file:

This is just out of curiosity but is the latter designed to be a backup of the original .qgs file?


texts:
I was shown the simplify function of https://mapshaper.org/, which removes unnecessary vertices in polygons. Now I'm looking for a similar function in QGIS and found GRASS generalisation. Since I'm not profound on generalisation algorithms and there are alot of them in GRASS generalisation could anybody tell me which algortihm I have to use to get an equal result?

What Mapshapers simplify does:

After you loaded a file you have a slider to simplify the polygon(s) The great part about is that it removes redundant vertices first and only after that starts to change the geometry. Additionaly the shared borders between polygons stay the same and you get no holes between them.


texts: Delineating catchments upstream of flow sites

texts:
I have many flow sites and I used them to delineate catchments upstream of each flow site using Spatial Analyst Tools > Hydrology > Watershed.

As you know using this tool, you'll get the catchments delineated as shown below (Figure 1).

Is there a way to get the catchments delineated as illustrated in Figure 2? In other words, as you move downstream, the catchment of each flow site will include all the catchments upstream of this site.


texts: Why do symbol categories disappear when unchecked in QGIS?

texts:
I am using Qgis 2.18.25

In a vector layer, when using the categorized style, when I uncheck a category the category disappears.

Is there any way to avoid this behavior?

Actually I do not really understand why this change was made.

Why would we want an unchecked category to disappear?


texts:
I would like to create circles around a point density layer that reflect the magnitude of the point density. I have created a point density estimation for a few towns. Each town varies in population, although the area does not change among them. I would like to create an "area of influence" around in town that reflects the population density. (Larger circles around higher populated towns.)

Is there a way to do this? I am using ArcGIS 10.1


texts: Unable to open password-protected raster in Georeferencer

texts:
I am using QGIS 3.4. I am unable to open a PDF file, a password-protected topo sheet in the Georeferencer. I do have the password.

The error message is:

  A password is needed. You can specify it through PDF_USER_PWD configuration option/ USER_PWD open option.

Where can I find this option?


texts: How to unckeck automatically the "No Data Value" predefined on raster background on QGIS 2.14?

texts:
I'm trying to import a bunch of rasters, but there is a parameter predefined on the TIFFs, where says that No Data Value should be marked. But I need this unchecked.

I opened the metadata, and found this lines:

<NoDataValue Sync="TRUE">256
</NoDataValue>

And even if I change that to FALSE, or remove completely the tag, or simply delete the files, the Rasters still are showing automatically with "No Data Value" checked.

This is the result I expect when I import the rasters:

But instead, all the rasters are been loaded automatically like this:

And I need to uncheck manually the box of "No data value" of every single raster like this:

What can I do to import all these rasters with no parameters of No Data Value?


texts: Retrieving User Input using ArcGIS API for JavaScript?

texts:
I am using the ArcGIS API for JavaScript to build a web application. Within the web app a geoprocessing service is executed which requires an input that is user dependent.

I thought of a text box where the user can enter a value, which should then be read by the application and used for the geoprocessing service.

How can this be done?
I searched on the ArcGIS API website and in the Internet, but I could not find any solution so far.


texts: Removing extra white space when exporting map from ArcMap

texts:
My map elements (map, legend, scale bar, and etc) are within the "Printer Margin", which I have set in the Page and Setup window.

However, when I export my map as a pdf, there is additional white space that is exported. The extra white space is only found at the bottom of the page not the top and the width of the exported map is fine.

How do I ensure that only the map elements within the Printer Margin are only exported without the extra whitespace?


texts: How to calculate the area of the convex hull of a points layer

texts: I have a shapefile in QGIS that I have created by importing a CSV of points, using Points2One to convert to a shapefile and then importing the shape. I am now trying to calculate the are of the shape but when trying to use the field calculator I can't figure out how to get the area? Can anyone help?


texts: Is there a way to identify features on OpenStreetMap.org?

texts:
Is there an identify features option for OpenStreetMap.org (OSM.org)?

Stemming from this question: How to get shapefile of river from OpenStreetMap? I'm looking for the ability to use OSM.org to figure out the tags and values for features I want to download through the Overpass API, or even QGIS for that matter.

If there was a way to use OSM.org to browse the data of OSM through the map interface, then click on a feature to get the information about it.

Currently, if you know what a feature is called, you can search for it, and you get the information I'm referring to as the output:

This would be the exact information I would expect if I were to use an Identify button and click on a feature on OSM.org.

This would greatly help, as the first referenced post suggests, the ability to pull down specific data from Overpass API or QGIS.

Thanks!


texts: DB manager in QGIS 2.0 UnicodeDecodeError

texts:
I followed this post:pgRouting 2.0 for Windows quick guide, and I want to use DB manager to load the test data into my PostgreSQL database.

I used url to get the test data:

I have set up the connection to my PostgreSQL 9.2, like this:

and the connection is fine, but as I fire up my DB manager, trying to connect to PostGIS, I got those errors:(errors got as I double clicked "PostgreSQL 9.2")

執行Python程式碼時發生錯誤:

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_model.py", line 404, in rowCount
    self._refreshIndex( parent, True )
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_model.py", line 456, in _refreshIndex
    if item.populate():
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_model.py", line 151, in populate
    if not connection.connect():
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\postgis\plugin.py", line 100, in connect
    return self.connectToUri(uri)
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\plugin.py", line 93, in connectToUri
    self.db = self.databasesFactory( self, uri )
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\postgis\plugin.py", line 67, in databasesFactory
    return PGDatabase(connection, uri)
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\postgis\plugin.py", line 134, in init
    Database.init(self, connection, uri)
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\plugin.py", line 172, in init
    self.connector = self.connectorsFactory( uri )
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\postgis\plugin.py", line 137, in connectorsFactory
    return PostGisDBConnector(uri)
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\postgis\connector.py", line 57, in init
    raise ConnectionError(e)
  File "C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\db_manager\db_plugins\plugin.py", line 40, in init
    msg = unicode( msg, 'utf-8' )
UnicodeDecodeError: 'utf8' codec can't decode byte 0xc4 in position 0: invalid continuation byte

Python 版本:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS 版本:
2.0.1-Dufour Dufour, d94c044

Python 路徑: ['C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins\processing', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python', 'C:\Users\Heinz/.qgis2/python', 'C:\Users\Heinz/.qgis2/python/plugins', 'C:/PROGRA~1/QGISDU~1/apps/qgis/./python/plugins', 'C:\PROGRA~1\QGISDU~1\bin\python27.zip', 'C:\PROGRA~1\QGISDU~1\apps\Python27\DLLs', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\plat-win', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\lib-tk', 'C:\PROGRA~1\QGISDU~1\bin', 'C:\PROGRA~1\QGISDU~1\apps\Python27', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg', 'C:\PROGRA~1\QGISDU~1\apps\qgis\python\plugins\fTools\tools']

Did I make some mistake in my process?

I am working under Windows 8.1, QGIS 2.0, PostgreSQL 9.2


texts: Debugging AttributeError: 'Result' object has no attribute 'save' from ArcPy script?

texts:
I Keep getting this error message when i run the code below, can anyone tell me why?

Error message:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File "C:\GEP_662\Final_Project\Final_Project1.py", line 65, in <module>
    result3.save('C:\GEP_662\\Final_Project\\satersym')
AttributeError: 'Result' object has no attribute 'save'

Code:

import arcpy

rasterLayer = "C:\GEP_662\\Final_Project\\TOC_IDW1.lyr"
rasterLayerName = 'TOC_IDW_LYR'

# Apply the symbology from the symbology layer to the input layer
outputraster = arcpy.ApplySymbologyFromLayer_management(rasterLayerName, rasterLayer)
outputraster.save('C:\GEP_662\\Final_Project\\rastersym')

print "symbology update done"


texts: pgr_dijkstra gives wacky routes sometimes with undirected graph

texts:
I am using the function below to find the shortest path between several origin-destination pairs, and it seems to work correctly most of the time.

CREATE OR REPLACE FUNCTION public.sp_od(
    orig integer,
    dest integer)
    RETURNS TABLE(shortest_path geometry)
LANGUAGE 'sql'

AS $BODY$

SELECT st_makeline(geom) as shortest_path
FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom, true) as cost FROM public."WA_roads"',
(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(ocentx as double precision), CAST(ocenty as double precision)), 4326) from all_trips_non_zero where origin = orig LIMIT 1) ASC
LIMIT 1),
(SELECT source FROM public."WA_roads"
 ORDER BY ST_StartPoint(geom) <->
     (select ST_SetSRID(ST_MakePoint(CAST(dcentx as double precision), CAST(dcenty as double precision)), 4326) from all_trips_non_zero where destination = dest LIMIT 1) ASC
LIMIT 1), directed := false
) as pt
JOIN public."WA_roads" rd ON pt.edge = rd.id;

$BODY$;

The above function takes an origin and destination zip code and using the corresponding lat, long for the zip code from 'all_trips_non_zero' table, finds the nearest "source" node in our road network closest to the origin and destination points to use for shortest path calculation.

One of the problematic path is shown in the image below.  Origin in red and destination is yellow. Brown is the network and blue is the shortest path from the above function.

For reference, the shortest path predicted by Google is as under.

I cannot explain some of the straight lines, not on the network, in the pgr_dijkstra path above. These are somehow part of the geometry and I get a total path length of around 360 miles, which is more than twice the distance predicted by Google. When I use directed:= true, I get no path for this OD pair. The "WA_roads" shapefile and "all_trips_non_zero" file are here. I further transformed the shapefile to EPSG:4326 as described here. The pgr_analyzegraph output is below:

NOTICE:  pgr_analyzeGraph('WA_roads',1e-06,'geom','id','source','target','true')
NOTICE:  Performing checks, please wait ...
NOTICE:  Analyzing for dead ends. Please wait...
NOTICE:  Analyzing for gaps. Please wait...
NOTICE:  Analyzing for isolated edges. Please wait...
NOTICE:  Analyzing for ring geometries. Please wait...
NOTICE:  Analyzing for intersections. Please wait...
NOTICE:              ANALYSIS RESULTS FOR SELECTED EDGES:
NOTICE:                    Isolated segments: 0
NOTICE:                            Dead ends: 214
NOTICE:  Potential gaps found near dead ends: 0
NOTICE:               Intersections detected: 1
NOTICE:                      Ring geometries: 2

Successfully run. Total query runtime: 1 secs 204 msec.
1 rows affected.


texts: "getting Unknown projection string error" serving EPSG:4326 tiles with open street maps mapnik renderd

texts:
The goal:
I am attempting to generate map tiles with EPSG:4326 instead of web mercator by specifying EPSG:4326 in the mapnik style file srs attributes

My current setup:
I am running Ubuntu 14.14 LTS and have followed the switch2osm.org manually-building-a-tile-server-14-04 instructions. Using the default projection, I was able to successfully import data and serve tiles to a client.  I am using the stack of osm2pgsql, renderd, mapnik, mod_tile, apache, postgres, and the OSM Bright stylesheet.
To import I used:
osm2pgsql --slim -d gis -C 16000 --number-processes 3 mapdata.osm.pbf

Then to run renderd:
sudo -u mapnik renderd -f -c /usr/local/etc/renderd.conf

At this point, I had everything working to serve up web mercator.

What I have completed to attempt to serve tiles in EPSG:4326:
First I reimported the data with osm2pgsql using the -l flag:
osm2pgsql  -l --slim -d gis -C 16000 --number-processes 3 mapdata.osm.pbf

I then updated the mapnik stylesheet map and layer srs attributes to match what was defined in http://spatialreference.org/ref/epsg/4326/mapnik/
srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

When starting renderd I will then receive the error:
"Unknown projection string, using web mercator as never the less.
 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

I also attempted a few other options:
srs="+init=epsg:4326"

srs="+proj=longlat +datum=WGS84 +no_defs"

After seeing a post recommending to check if proj-epsg package is installed, I was heading down the path that proj was missing on my system  (https://github.com/mapnik/mapnik/wiki/InstallationTroubleshooting#proj_init_error)
Per this proj init error guide I reinstalled proj to no effect.
I confirmed I do have it:
~$ proj
Rel. 4.8.0, 6 March 2012
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

I have checked my /usr/share/proj/epsg file and it contains an entry:
# WGS 84
<4326> +proj=longlat +datum=WGS84 +no_defs  <>

The guide shows a different path to /usr/local/share/proj/epsg.  On my system it is /usr/share/proj/epsg

Question:
Is srs="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" to output EPSG:4326 a valid option to generate tiles with renderd or is this not a possibility?

For reference, current configuration files:
renderd.conf:
[renderd]
socketname=/var/run/renderd/renderd.sock
num_threads=4
tile_dir=/var/lib/mod_tile
stats_file=/var/run/renderd/renderd.stats

[mapnik]
plugins_dir=/usr/local/lib/mapnik/input
font_dir=/usr/share/fonts/truetype/ttf-dejavu
font_dir_recurse=1

[default]
URI=/osm_tiles/
TILEDIR=/var/lib/mod_tile
XML=/usr/local/share/maps/style/OSMBright/OSMBright.xml
HOST=localhost
TILESIZE=256

problem snippet EPSG:4326 version of mapnik OSMBright.xml:
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map[]>
<Map srs="+proj=longlat +datum=WGS84 +no_defs"
font-directory="./fonts" background-color="#c4dff6"
maximum-extent="-20037508.34,-20037508.34,20037508.34,20037508.34">

working snippet web mercator version of mapnik OSMBright.xml:
<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map[]>
<Map srs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs +over"
font-directory="./fonts" background-color="#c4dff6"
maximum-extent="-20037508.34,-20037508.34,20037508.34,20037508.34">


texts: Understanding Number of Feature Classes vs Connection Time with Enterprise Geodatabases

texts:
So here's the situation - at the moment I am having a problem where the more feature classes I load into a ArcGIS 10.2 enterprise geodatabase in PostgreSQL, the slower the connection time becomes.  I have gone from 40 seconds with 2000 feature classes (both vector and raster) to 2 minutes with 6000 feature classes.

Has anybody else noticed this problem when adding more feature classes, if so did you manage to overcome the problem, and how did you fix it?


texts:
I am looking for an up-to-date communes administrative boundraries (LAU) of several countries in Europe (France + adjacent countries). I need it to be quite precise (less than 5m), since I intend to use it to retrieve commune names from a point (bus stop). I was thinking about OpenStreetMap and Overpass Turbo, but the coverage seems to be too wide to be processed, even country by country, here is one of my resquest :

[out:xml][timeout:600]; (
  area['admin_level'='2']['name'='France'];
  rel['boundary'='administrative']['admin_level'='8'](area);
);
(._;>;);
out;

I also found

[something interesting] on the eurostat website, but it is from 2013 ...
[Official LAU], but it not precise enough for my needs
[GDAM data], but again, not precise enough


texts:
I am exporting a grass vectior dataset to a Postgis dump file using v.out.ogr (reformatted for clarity)

v.out.ogr -a -n input=my_layer@mapset\
                output=my_sql_file.sql format=PGDUMP

It exports using insert statements by default, but the ogr2ogr docs suggest that using a switch enables COPY format instead (--config PG_USE_COPY YES).  COPY is much faster and is also easier to debug when errors occur.

I tried using PG_USE_COPY=YES as both a dsco and lco option, and exporting it (ubuntu 18.04 command line in grass environment).

export PG_USE_COPY=YES

Is it possible to export using COPY format?

What am I doing wrong?

[edit]

My problem was caused by creating a layer by merging server datasets, one of which had an ogc_fid field (default primary key for ogr2ogr).  This caused loading to fail after several thousand insert statements.  I solved it by editing the sql file to rename the primary key and add ogc_fid as an integer field.


texts: Calculate Locations not adding Status field

texts:
When I run the Calculate Locations tool, it inconsistently adds the Status field.  It consistently adds the SourceID, SourceOID, PosAlong, SideofEdge, SnapX, SnapY, and Distance fields.

I've tried the tool from both Toolbox and Model Builder, and both tools will apply the field to the same point feature class input (e.g. fire stations), but neither will apply it to another feature class input (e.g. addresses).

Cross-posted for visibility at GeoNet


texts: How to create multiple buffers around line at equal intervals without typing in every value?

texts:
I have a line in ArcMap 10.2, and want to create lines to the right of it at 100 m, 200 m, 300 m etc. up until 1.7 km. Obviously typing in these values under 'multiple ring buffers' will take too much time. Is there a way for ArcMap to compute this automatically across the equal intervals?

If you could answer through basic ArcMap methods that would be great. I do not understand the coding language that some of you use!


texts: MMQGIS converts very large text numbers to a numerical missing value

texts:
I used MMQGIS plugin to change a table of numbers in text format (UN world economic data) which had large values of items like GDP 132000000000 (132 billion) were all converted to -2147483648 which I suppose is the numerical value for missing value or similar. 132 million converts OK but not 132 billion.

The image above illustrates the problem. I started with the table on the right; table on the left is the reslt. I am using QGIS 1.8.  How can I fix this?


texts: Can QGIS get the GPS Coordinates via the 'Windows Location Provider'?

texts:
Can Quantum GIS get the GPS Coordinates via the 'Windows Location Provider' in Windows 7/8 instead of reading a virtual com port?

If not, could this be developed as a Windows-specific plugin?  (Is there a precedence for developing Windows-specific plugins?)


texts: How do you calculate the area of polygons which overlap a buffer around a point?

texts: I'm new to QGIS (Dufour 2.0.1) and have plotted nests as points on a googlemaps base-map and have created buffer zones around the nests. I've also plotted ponds as polygons on the same base-map. I want to know how to calculate the area of the ponds which fall within the buffer zone PER NEST.


texts: How to enable projection transformation in QGIS

texts:
My QGIS project has two layers, one layer is retrieved by Google layer plugin which fetches the static Google map image whose CRS is +proj=merc +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs and another vector layer is a postgis layer whose CRS is +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs.

I digitalized some polygon features from Google Layers and added it to the vector layer. However, QGIS seems do not make the automatic CRS transformation. How can I change its projection?

Thanks


texts:
Given that an area of interest is partially cloud covered, how can you mosaic two separate collections to get a single cloud free mosaic?

Given I have a S2 image over a study area, that has a cloudy eastern half:

var collection = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2017-01-01', '2017-03-01')
.filterBounds(geometry5)

and imagery over the same region from a previous date that has a cloudy western half:

var collection2 = ee.ImageCollection('COPERNICUS/S2')
.filterDate('2016-11-13', '2016-12-23')
.filterBounds(geometry4)

Is there a way to mosaic the data together to get a single collection/dataset?

This is a straightforward task with gdal but I can't find it in the Earth Engine documentation.

As a visual example:

I have this collection:

and this separate collection:

and I want to merge/mosaic them into a single collection:


texts: Permanently attaching Image file to Feature Class in 10.1

texts:
I need to forward on a File Geodatabase, and I am using the attachment manager to attach the file images associated with the different features in my Feature class.

If I forward on this FGD, will the images be permanently there?  I have not been able to find anything about whether it is permanent, or I would need to send the source files along, and the attachments would have to be redirected on the end of the recipient?

If not is there another way to permanently attach an image file to a Feature class in a File Geodatabse?  These photos are not geotagged, nor do they have a time stamp on them so neither of the photo tools would work on them.  There is no possibility that this can be in SDE.

Thanks.


texts:
I used the following code:

this.installDrag({layer:this.poiLayer,opt:'load',id:'poiLatlong'});
this.installDrag({layer:this.addNewPoiLayer,opt:'add',id:'poiLatlong'});

but drag is installing for only for the last layer I call. How can I apply installDrag to more than one layer?


texts: How to determine if a GeoTIFF image is a photo or a raster map using GDAL and Python

texts: Does anyone have a method to determine if a GeoTIFF image is a photo (e.g aerial ortho) vs a raster map (e.g topographic map produced from vectors) using GDAL and Python? Both images are RGB/RGBA uncompressed with a byte datatype.


texts: Is there a way to scale to 1:1 in QGIS Cloud?

texts:
Someone knows if with QGIS cloud Pro Plan, it's possible have a scale till 1:1. I use QGIS to draw small stuff on my gis (archaeological finds, as really small pottery's sherds) and I need a scale smaller than 1:624 as now I can view with QGIS cloud free plan.

Thank a lot and best regards


texts:
I've been dealing with various problems after re-installing QGIS and other GIS-related software (see QGIS corrupted plugins upon reinstallation).

The solution worked for a bout a day and now I lost the processing tool again. An other error message appeared when opening QGIS:

Impossible to load processing plugin provokes error when calling classFactory()

UnicodeDecodeError: 'ascii' codec can't decode byte 0xe9 in position 18: ordinal not in range(128)

Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 306, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/F�lix/.qgis2/python/plugins\processing\__init__.py", line 37, in classFactory
    from processing.ProcessingPlugin import ProcessingPlugin
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/F�lix/.qgis2/python/plugins\processing\ProcessingPlugin.py", line 47, in
    if cmd_folder not in sys.path:
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 58, in showWarning
    "warning:%s\ntraceback:%s" % (warnings.formatwarning(message, category, filename, lineno), stk),
UnicodeDecodeError: 'ascii' codec can't decode byte 0xe9 in position 18: ordinal not in range(128)

Also, I have a second error when trying to re-install the processing extension through QGIS.

ImportError: cannot import name defaultOutputFolder

Traceback (most recent call last):
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 281, in loadPlugin
    __import__(packageName)
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/F�lix/.qgis2/python/plugins\processing\__init__.py", line 28, in
    from processing.tools.dataobjects import *
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/OSGEO4~1/apps/qgis/./python/plugins\processing\tools\dataobjects.py", line 35, in
    from processing.core.ProcessingConfig import ProcessingConfig
  File "C:/OSGEO4~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
  File "C:/Users/F�lix/.qgis2/python/plugins\processing\core\ProcessingConfig.py", line 32, in
    from processing.tools.system import defaultOutputFolder
ImportError: cannot import name defaultOutputFolder

Any ideas on how to solve this one?


texts:
I am trying to get a centroid for several polygons on a global scale, in a way that it would correctly represent the common centre of all the polygons, taking Earth's curvature into account. The polygons are in the WGS 84 coordinate system.

For testing purposes, I created a test shapefile with two polygons: one that encompasses Greenland, and one that encompasses the area east of Mongolia up to the sea. These areas are quite a bit to the north, and their real common centroid is supposed to be somewhere close to the North Pole. Probably around the Laptev Sea, since Greenland is a bit smaller area-wise (although you couldn't tell that from the default projection QGIS is using to display WGS 84).

If I use the "Polygon centroids" built-in algorithm, the result is somewhere around St. Petersburg. That is way too south-west from where it's supposed to be.

Another test created even more hilarious results. The two areas I selected that time were Alaska and the Kamchatka peninsula. The centroid should be at the very north of the Pacific ocean. But the result of "Polygon centroids" is actually in Nunavut, Canada, which is on the opposite side of where it should be. It's put there because Alaska is bigger than Kamchatka, otherwise it would be put somewhere in the Atlantic ocean.

One way of getting more correct centroids I could see is by changing the projection system, for instance, to an Azimutal Equidistant one, where the reference is somewhere close to the North Pole. However, if I do change the projection system, the polygons are still in WGS 84, and QGIS still outputs the same incorrect results if I do another pass of "Polygon centroids". And if I set the polygon layer's CRS to the projected system, the old coordinates are retained instead of being recalculated for the new projection.

This is how both test cases look like when projected, and shows clearly how incorrect the results are (the centre of the image is the North Pole):

So perhaps there is a way to permanently reproject the polygons to the projected CRS? Or some other way of getting correct polygon centroids?


texts: How do I convert my coordinates to the limit of geographic region (outer edge)?

texts:
I am trying to use r.in.bin. The parameters of north, south, east, and west needs values:
north=float
Northern limit of geographic region (outer edge)
south=float
Southern limit of geographic region (outer edge)
east=float
Eastern limit of geographic region (outer edge)
west=float
Western limit of geographic region (outer edge)

I only have the values for my map, which is in this format:
north: 15N
south: 10N
east: 125E
west: 120E
rows: 6000
cols: 6000

How do I convert them?


texts: Point features appear in PDF export when they shouldn't

texts:
I'm using QGIS and have some issues when exporting maps using point layers.

In this case, I've simply displayed proportional circles from a point layer using the size assistant, setting their minimum size to 0. The circles are set with a thin outline.

The result looks fine on the project itself and on the PNG export: the points whose value equals 0 can't be seen. Thus there's no issue with their outline I assume.

However, these points, which shouldn't be displayed, are still visible on the PDF export.

The same issue can occur when point layers are solely used to display labels (with the label property), in which case the underlying point should be invisible.

Any idea on how to fix this?


texts:
After adding the following to the top of my index.html page

    <!DOCTYPE html>

I begin to receive the following error on my GetMap WMS request

    Resource interpreted as Image but transferred with MIME type application/vnd.ogc.se_xml

Any thoughts on why this would occur are appreciated.  Here are the URL params of my WMS request:

    BBOX=-15130249.425402,2538320.1060898,-5869750.5745982,6461679.8939102
    FORMAT=image/png
    HEIGHT=802
    LAYERS=EFETAC-NASA_1YrBaseline_previous1
    MAXEXTENT=NaN,,,
    PROJECTION=EPSG:900913
    REQUEST=GetMap
    SERVICE=WMS
    SRS=EPSG:900913
    STYLES=
    TRANSPARENT=true
    UNITS=m
    VERSION=1.1.1
    WIDTH=1893


texts: Defining coordinate system for every GeoDataFrame (shapefile) in list

texts:
I am looking to set a coordinate system for every GeoDataFrame in a list of them. So I need this:
gdfs[0].crs = "EPSG:6668"
gdfs[1].crs = "EPSG:6668"

And so on but iteratively. I tried:
([gdfs.crs = "EPSG:6668"] for file in gdfs)

And:
for file in gdfs:
    gdfs.crs = "EPSG:6668"

But not working, any solution?


texts:
First, I'm trying to call a function from a module (final_functions) I created in a separate script. However I receive the following error:

Traceback (most recent call last):   File "C:\Python27\ArcGIS10.3\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 325, in RunScript
    exec codeObject in __main__.__dict__   File "C:\EsriPress\Python\Data\Final\Test.py", line 8, in <module>
    check = final_functions.FieldExists(fcname, fieldname)   File "C:\EsriPress\Python\Data\Final\final_functions.py", line 6, in FieldExists
    namelist = [] NameError: global name 'fc' is not defined

Here's the main script:

import arcpy
from arcpy import env
arcpy.env.workspace = "C:/EsriPress/Python/Data/Final/Scratch.gdb"
import final_functions
fclist = arcpy.ListFeatureClasses('Fire*')
for fc in fclist:
    fc = fcname
    fieldname = "SFCount"
    check = final_functions.FieldExists(fcname, fieldname)
    print check

Here is the module (final_functions.py) I'm importing:

def FieldExists(fcname, fieldname):
    namelist = []
    fields = arcpy.ListFields(fcname)
    for field in fields:
        namelist.append(field.name)

    if fieldname in namelist:
        return "Field: SFCount already present in Feature Class:     {0}".format(fcname)
    else:
        arcpy.AddField_management(fcname, fieldname, "TEXT")

if __name__ == "__main__":
    fclist = arcpy.ListFeatureClasses('Fire*')
    for fc in fclist:
        fieldname = "SFCount"
        fcname = fc
        checkfields = FieldExists(fc, fieldname)
        print checkfields

When I run final_functions on its own, it works great. I'm only given the error when trying to call it in the main script.

And if you're wondering, I'm just checking a list of feature classes for the existence of a field, and adding a new field if it doesn't exist.


texts: Known Intersecting Polygons returning false for .intersects() in geopandas

texts:
I have 2 geodataframes:

An original shapefile of water features read in as a geodataframe. (playas)
A gdf as a result of a convex_hull operation with only polygon geometries. (polygons)

There should be intersection between these two layers, but the following code returns an empty geodataframe:

aff_playas = playas[playas.intersects(polygons)]

I am trying to think of what the problem could be with these two layers. I believe it is one of the following:

A projection inconsistency with the created hulls_gdf
function failure due to how detailed, varied the playas layer is.

The two layers plotted below have similar extent and are verified to intersect in Qgis, ArcGIS. The third is the output from the intersects().


texts: Can't install ArcGIS License Manager 10.1 on Win Server 2003

texts:
i try to install a License Manager on a Windows 2003 Server.

The setup starts, everything looks fine. When i press "next" the setup crashes, the setup window closes. No error. No entry in the event viewer...

I already tried to create a installation log.

  "C:\MyPackage\Setup.exe" /L*V "example.log"

But nothing is created...

Whats the problem? :/


texts: Is there a way to independantly hillshade the overviews of a GeoTIFF pyramid in GDAL?

texts:
Is there a way to get gdaladdo, gdaldem or another GDAL tool to create a  hillshade where each overview is shaded from the corresponding overview in the DEM rather than by resampling a single full size hillshade that's computed once?

A way to build turn a collection of different resolution images of the same extent into a single GeoTIFF with overviews would be sufficient.


texts: Walking route that covers the largest area in shortest distance

texts:
I regularly need to walk across field searching for bird nests. The field are often squares or polygons and the bird nests could be anywhere on the ground in the field. To stand the best chance of finding bird nests, I generally zig-zag across the fields. I have also assumed that zig-zagging covers the largest area in the shortest distance.

Below is a hypothetical field and my current zig-zag approach for walking the field. However, I'm wondering if there is a better way of walking the field in order to cover the largest area?

Assume that I can see 10m either side of line. Can anyone provide R code for routes that would cover 10%, 20%, 30%……90%, 100% of the field?

I'm particularly interested in R code to calculate a walking route, hence why I've used the R tag.


texts:
I have a land use raster file and am trying to use the raster calculator to query the file for surface water.

The saga formula (11 = surface water band):

(a, 11)

This does not work.

Can anyone help me with the raster calc syntax to perform this?


texts: How to change the values of the push buttons in QT designer

texts:
I am using the standard dialog box with ok and cancel options.i Have also added one more button exit.

self.dlg.show()
result=self.dlg.exec_()
print 'The result is',result

it is giving the result of 1 for ok. but it is giving the result of 0 for both cancel and the exit buttons. I want to assign different functionalities to both the buttons. So how to chnage the result values for both the buttons?
like i want to print the result value of exit to be 2.


texts: Filtering footprints to given bounding box

texts:
I'm trying to get footprints with landuse tag within sample bounding box with query below:

[out:json][timeout:180];
(
  way(poly:"50.859515 18.568600 50.859515 18.580753 50.862946 18.580753 50.862946 18.568600 50.859515 18.568600")["landuse"];(._;>;);
  relation(poly:"50.859515 18.568600 50.859515 18.580753 50.862946 18.580753 50.862946 18.568600 50.859515 18.568600")["landuse"];(._;>;);
);
out body;

Unfortunately, returned ways are out of my bounding box, as shown on the picture below (bbox is approximately red square) so the question is how to filter results to retrieve data only in my bounding box range?


texts: What transformer should I use for `gdal.Warp`?

texts:
The documentation for gdal.Warp (http://gdal.org/python/osgeo.gdal-module.html#Warp) lists the following alluring options:
  tps --- whether to use Thin Plate Spline GCP transformer
  rpc --- whether to use RPC transformer
  geoloc --- whether to use GeoLocation array transformer
  polynomialOrder --- order of polynomial GCP interpolation

I'm having trouble finding information about these options. Is there a reason one would want to choose one over the other?


texts: How to merge two RGB raster tables with overlap by max in each band in postgis?

texts:
I need to merge RGB raster tables from Sentinel-2 data into a mosaic of same date using ST_Union(), but following the example of Postgis manual for all bands and all input rasters like this one I got errors.

Ps.: All rasters were processed to ressample and reproject to same resolution and projection before to upload at postgresql server.

--My code

select ST_Union(ARRAY[rast1, rast2], 'MAX') from
    sentinel.rasterTable1 as rast1,
          sentinel.rasterTable2 as rast2;

-- the error
ERROR:  function st_union(sentinel.rasterTable1, sentinel.rasterTable2) does not exist

How can I fix this to get the max by each band?

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

First appoach by @Pierre using this code:

 create table mergedData as SELECT ST_Union(rast, 'MAX') rast FROM
 (SELECT rast FROM rasterTable1
  UNION ALL
  SELECT rast FROM rasterTable2) foo

Aparently, only copy the first image listted in the code above.


texts: Using wildcards to work around GDAL character limit

texts:
I am attempting to use the 'Merge' & 'Create VRT' tools in QGIS 2.0 on
a folder with 3000+ raster tiles in.

However whenever I try to perform this operation the following error appears:

  The process failed to start. Either the invoked program is missing, or
  you may have insufficient permissions to invoke the program.

The process works when I select only a few tiles, so I can only assume it's due to a character restriction due to the number of files?

Upon searching the website for a solution I came across the concept of using 'Wildcards' to get around this problem:

Having problems creating VRT in QGIS from png files

However I cannot find any clear explanation on how to use them, and I don't have much experience with using OSGeo4W outside of the QGIS GUI.


texts: Slovenian Municipalities - right spelling

texts: I'm working on Slovenian Municipalities. Matching them with georeferencing only hits around 80 pct. I have then tried MunicipalityId, but it doesn't work at all. Where can I get the list of names on municipalities, used by CartoDB, so I can make a reference dataset?


texts: Style MapServer layers easily (or as easy as possible)

texts:
I am working with MapServer WMS in the form of mapfiles and consuming them into a OpenLayers map.

I am trying to find out a way to ease the symbolization process of the various layers in a mapfile.

I have not started diving into the STYLE directive yet, but before I do that I'd like to know whether there's an easy way to (let's say) preview the style of my layers (an application?) and export it into a (e.g.) SLD which the MapServer can read.

I don't know if this is the best (or feasible) way to go about styling mapfile's layers, but I am open to any kind of suggestion.


texts: I am looking for a document which explains what the API of QGIS is, with examples of this API but in French.


texts: SQL procedure(s) to convert GPS traces to an ordered set of junctions in a network

texts:
I have three datasets that I have loaded into a geodatabase (spatialite or PostGIS):

A set of GPS points collected from individual trips through urban areas. Points have a unique point id, a group trip id, as well as the usual GPS attributes (speed, accuracy, elevation, timestamp, etc).
A set of points corresponding to junctions in a road network.
A set of polylines corresponding to the road segments connecting junctions. These can be identified as a tuple of junction IDs.

I am trying to relate the GPS points back to routes in a network. Ultimately, from each set of GPS points constituting a trip, I would like to identify an ordered set of junctions from my junctions dataset. (Order, JunctionID)

Initially, I started by creating a modestly sized buffer around my junctions, and then simply intersecting the GPS points, grouping by trip id. This gets surprisingly close to identifying a reasonable set of junctions per trip. However it has a few deficiencies:

Where junctions are close together, individual GPS points can be captured multiple times. Some tie-breaking mechanism is needed.
Similarly, where the urban canyon effect comes into play, there may be stray points erroneously by other junctions.
There is the possibility that an individual backtracked on their route, meaning that a junction could have been visited multiple times.

I am envisioning various constraints that could potentially be used to minimize junction selection error as well as to identify cases of backtracking. Those constraints could involve choosing junctions with the most GPS points, or only pairs of junctions that have a segment connecting them. I could also imagine looking at the timestamps to identify instances of backtracking.

Has anyone done something similar to this who can spare me reinventing the wheel? Or if I'm completely barking up the wrong tree, and there is a Better Way of identifying a network route from GPS traces in a database-driven way, I'd be open to that too!


texts: Integrating Landsat ETM+ and OLI surface reflectance data

texts:
I want to do a time-series analysis on Landsat surface reflectance data for an extended period of time (2005-2018), which requires integration of Landsat 5, 7, and 8. Is there any procedure to make the data from ETM+ and OLI sensors compatible so they can be used in a single time series together? I looked at the scholarly publications but only found one paper dated 2014, suggesting a linear regression technique for converting OLI to ETM+.

https://www.researchgate.net/publication/277674991_Continuity_of_Reflectance_Data_between_Landsat-7_ETM_and_Landsat-8_OLI_for_Both_Top-of-Atmosphere_and_Surface_Reflectance_A_Study_in_the_Australian_Landscape


texts:
Given a database of GPS points, a center GPS location, and a radius R (in meters), how can I find all points in the database that are within a distance of R of the center point?
Will the distance iso-lines be an ellipse?


texts: Getting last image dates automatically in Google Earth

texts:
I want to check whether Google Earth has an up-to-date image (taken in this year) for a given location. I can do it manually by checking the location on Google Earth Pro. However, since I am interested in many location at the same time, I am looking for a table that shows cities, for example, and their corresponding last image date. With this information, I want simply sort cities that has an image in 2020, for example, and analyze these cities.
Is there any way of getting it?


texts: how to divide a shapefile and its data into tiles?

texts:
I have a population map of counties shapefile. I need to divide this data into tiles in order to have a gridded population map.

population_of_cell = population_of_county / area_of_county * area_of_cell

in grids intersecting two or more polygons, respecting the area share of each polygon, the population of cell is determined.
how can I divide so using QGIS?


texts: Setting elevation data for tie points in Pix4D

texts:
I am currently working in Pix4D and am trying to assign elevation data to my selected tie points. I have aerial images taken using a camera with GPS from an aircraft and am working with a fairly rectangular, level area of interest (AOI). I want to set a base elevation for the AOI, and have selected 6 sets of tie points - 1 at each corner and then 2 on either sides. I will then be giving them all the same elevation.

I have carried out the processes as detailed in the following post - Reset a base elevation in Pix4D? - but all I was given was the option to change the altitude of the what looks like the camera (which is in the aircraft).

Does anyone have any idea how to assign elevation data to tie points specifically?


texts:
I would like to know if any of you has shapefiles of Berlin?

I have been using https://fbinter.stadt-berlin.de/fb/index.jsp but it doesn't allow the layers in shapefile format and also https://opendata-esri-de.opendata.arcgis.com/search?sort=-modified but I can't find the shapefile of Groundwater.

I am working on my thesis project and need this information.


texts: Conditional calculation in QGIS Raster Calculator

texts:
Is it possible in QGIS either directly or via plugins to do a raster calculation including an if statement?

e.g., if raster1 > 0.3, then raster 1 * raster 2, else raster1 x raster3... something like that?


texts: bbox for a particular feature of a layer in geoserver

texts:
I am trying to get the bbox for a particular feature for a layer from my geoserver. For this I am trying using the following url:

http://example.com/geoserver/wfs?version=1.3.0&request=describeFeatureType&outputFormat=application/json&service=WFS&typeName=layerName&PROPERTYNAME=(DISTNAME,Sylhet)

Ok fine, I am getting the bbox in the response. But the response gives me the coordinates for that particular feature too. Its very clumsy to me and I just want the bbox as response not the coordinate.

For this how can I change my request??
Please help.


texts: Creating Tin from Elevation Points in QGIS?

texts:
I have a shapefile of points with values that represent elevation at the given point. I know that I can create a TIN from these points, but the fields in the TIN do not represent the elevation. Rather, they represent an ID field for each node in the triangle.

I can also use the Interpolation plugin with Triangulated interpolation as the method. However, in this case it produces a DEM.

How can I produce a TIN where each node in the TIN represents the elevation at that point?


texts: Proj4 WKT - Cartesian with origin to EPSG900913

texts:
I'm converting wgs84 to cartesian with origin using simple function in c#

public static double[] LL_TO_XY(double latitude, double longitude, double originLatitude, double originLongitude)
    {

            //var xx,yy,r,ct,st,angle;
        double xx;
        double yy;
        double r;
        double ct;
        double st;
        double angle;
        double rotation_angle_degs = 0;
        double xoffset_mtrs = 0;
        double yoffset_mtrs = 0;
        double pxpos_mtrs;
        double pypos_mtrs;
        double[] p = new double[2];

            angle = DEG_TO_RADIANS(rotation_angle_degs);

        xx = (longitude - originLongitude)*METERS_DEGLON(originLatitude);
        yy = (latitude - originLatitude)*METERS_DEGLAT(originLatitude);

                r = Mathd.Sqrt(xx*xx + yy*yy);

                /* alert('LL_TO_XY: xx=' + xx + ' yy=' + yy + ' r=' + r);
    return false;*/

                //if(r)
                //{
                    ct = xx/r;
                    st = yy/r;
                    xx = r * ( (ct * Mathd.Cos(angle)) + (st * Mathd.Sin(angle)) );
            yy = r * ( (st * Mathd.Cos(angle)) - (ct * Mathd.Sin(angle)) );
                //}
                pxpos_mtrs = xx + xoffset_mtrs;
                pypos_mtrs = yy + yoffset_mtrs;

        p [0] = pxpos_mtrs;
        p [1] = pypos_mtrs;
        return p;

    }

Now I want to convert these cartesian values into EPSG900913. I can use proj4 for this so I think it should be better.

I know that to convert from WGS84 to EPSG900913 I need to define WKT for proj4 like this:

private static string wktEPSG900913 =
    "PROJCS[\"WGS84 / Simple Mercator\", " +
    "GEOGCS[\"WGS 84\", " +
    "DATUM[\"World Geodetic System 1984\", SPHEROID[\"WGS 84\", 6378137.0, 298.257223563,AUTHORITY[\"EPSG\",\"7030\"]], " +
    "AUTHORITY[\"EPSG\",\"6326\"]]," +
    "PRIMEM[\"Greenwich\", 0.0, AUTHORITY[\"EPSG\",\"8901\"]], " +
    "UNIT[\"degree\",0.017453292519943295], " +
    "AXIS[\"Longitude\", EAST], AXIS[\"Latitude\", NORTH]," +
    "AUTHORITY[\"EPSG\",\"4326\"]], " +
    "PROJECTION[\"Mercator_1SP\"]," +
    "PARAMETER[\"semi_minor\", 6378137.0], " +
    "PARAMETER[\"latitude_of_origin\",0.0], " +
    "PARAMETER[\"central_meridian\", 0.0], " +
    "PARAMETER[\"scale_factor\",1.0], " +
    "PARAMETER[\"false_easting\", 0.0], " +
    "PARAMETER[\"false_northing\", 0.0]," +
    "UNIT[\"m\", 1.0], " +
    "AXIS[\"x\", EAST], AXIS[\"y\", NORTH]," +
    "AUTHORITY[\"EPSG\",\"900913\"]]";

What WKT I need to use to convert from cartesian with origin to EPSG900913?


texts: Different result featureCount() on QgsVectorLayer QGIS 1.8 vs QGIS 2.0?

texts:
I use the following code in a layer action:

<uri definition etc>
sql = 'id = "[% "id" %]"';
uri.setSql(sql);
layername = 'TableName [% "id" %]';
vlayer = QgsVectorLayer(uri.uri(),layername,'spatialite');
<some code> if (vlayer.isValid() and vlayer.featureCount() > 0) else <some code>

The SQL-query gives as result 0 or more features. vlayer.featureCount() returns 0 or more as expected in QGIS 1.8, but always 359 in QGIS 2.0 (total features in Spatialite table with name TableName), so my code breaks in QGIS 2.0.
When the SQL-query has no result and I add vlayer to the interface and turn Show Feature Count on the layers panel shows two counts (?!) [359] and [0]. When the query returns 1 feature, these counts are [359] and 1. Open Attribute Table shows 0 and 1 record which is all fine. So the SQL-query seems to run fine, but featureCount holds the wrong number in my opinion.

Any clue?


texts: How to use SpatiaLite functions in a Python script?

texts:
I'm trying to find the intersection coordinates between a line (knowing its start and end points) and a polygon (SQlite file).

I've managed to make it work in QGIS by creating a shapefile for the line and then using DBmanager and the following SQL statement:

SELECT astext(st_intersection(polygon.geom, line.geom)) as XY
from polygon, line
WHERE st_intersects(polygon.geom, line.geom)

It works just fine but now I'm trying to make it work in a python script with the XY coordinates of the line and the polygon.sqlite file as inputs.

I'm trying to write the code and I've come to this so far:

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

import gdal
import ogr
import sqlite3
import math
import shapefile
import pyspatialite

fichier=open("C:\\MyPath\\polygon.sqlite","r")
connR =sqlite3.connect("C:\\MyPath\\polygon.sqlite")

curR =connR.cursor()

#trying to add a table for the line:

curR.execute("""
               CREATE TABLE my_line(
                   id INTEGER PRIMARY KEY,
                   XMin INTEGER,
                   YMin INTEGER,
                   XMax INTEGER,
                   YMax INTEGER,
                   geom MULTILINESTRING
                )
""")
connR.commit

# user coordinates inputs

XA=input('X A = ');
YA=input('Y A = ');
XB=input('X B = ');
YB=input('Y B = ');

curR.execute("""INSERT INTO route VALUES (?,?,?,?,?,?)""", (1, XA, YA, XB, YB,"MULTILINESTRING"))
connR.commit()

#trying to use SpatiaLite functions:

curR.execute("""
        SELECT astext(st_intersection(polygon.geom, line.geom)) as XY
        from polygon, line
        WHERE st_intersects(polygon.geom, line.geom)
""")
connR.commit()

As you can see, my first problem is that I don't know how to create a "line" object in the table 'my_line'.
My second problem is that I don't know how to use the SpatiaLite functions in a python script.


texts:
I need to split a huge raster into smaller parts to be able to process them part per part. (It's not about the size, however the extent is spread across several latitude degrees which is a problem for solar irradiance calculations.) Therefore, I want to extract these parts from the Raster, using zonal Shapefiles.

Is there a tool in ArcGIS to create repeating rectangular shapes?


texts: Mapbox.js setStyle differently for each element

texts:
I'm using geoJson in a featureLayer so I can use the setFilter-Function. My question is: How can I use the setStyle-Function and give Elements with different properties different colors?

My Code:

var jsonLayer = L.mapbox.featureLayer(geoJson).addTo(map);
//Filter for start_year
jsonLayer.setFilter(function(f){return f.properties["start_date"] <=outyear;})
jsonLayer.setStyle(getStyle(this));

function getStyle(feature){return{
stroke: '#555555',
strokewidth: 1,
fillColor: getColor(feature)}};

function getColor(feature){
    if (feature.properties.start_date<="1700"){return "#c78432"};
    else if (feature.properties.start_date>="1800"){return "#F9A63F";}
    else {return "#B8B8B8"};}

This obviously does not work correctly.