texts:
I have this query

SELECT
data_type, data_time, model_name, ST_Value(raw_data, 1, 202,417) AS val_0
FROM data
WHERE data_type = 'type_a'
  AND data_time >= '2015-12-06 00:00:00'
  AND data_time <= '2015-12-07 00:00:00'
ORDER BY
data_type DESC, data_time ASC

It runs about 80ms (data time sampling is every 3 hours => result contains 8 rows). Can this query be speed up?

Here is output from explain:

"Sort  (cost=14.46..14.48 rows=8 width=39)"
"  Output: data_type, data_time, model_name, (st_value(raw_data, 1, 202, 417, true))"
"  Sort Key: data.data_time"
"  ->  Bitmap Heap Scan on public.data  (cost=4.87..14.34 rows=8 width=39)"
"        Output: data_type, data_time, model_name, st_value(raw_data, 1, 202, 417, true)"
"        Recheck Cond: ((data.data_time >= '2015-12-06 00:00:00'::timestamp without time zone) AND (data.data_time <= '2015-12-07 00:00:00'::timestamp without time zone) AND (data.data_type = 'type_a'::text))"
"        ->  Bitmap Index Scan on data_pk  (cost=0.00..4.87 rows=8 width=0)"
"              Index Cond: ((data.data_time >= '2015-12-06 00:00:00'::timestamp without time zone) AND (data.data_time <= '2015-12-07 00:00:00'::timestamp without time zone) AND (data.data_type = 'type_a'::text))"

Database structure:

CREATE TABLE data
(
  raw_data raster,
  data_time timestamp without time zone NOT NULL,
  data_type text NOT NULL,
  model_name text,
  CONSTRAINT data_pk PRIMARY KEY (data_time, data_type)
)

Edit:

Raster has 1 band and there is 1 row for every hour (so each row is one image). Raster resolution os 1500x750, top left corner is set to (0, -90), skew is 0 and scale 0.25 in both directions.


texts: Filter a date interval (month) shapefile

texts: I'm trying to filter some data of a shapefile (with "Filter") for QGIS to show points in a date interval of a shapefile. I have a shapefile which has a table of attributes of sightings of animals with time and date and other attributes relevant for me. I'm now trying to create shape files of each month, I do a query using the "Filter" and I can't seem to do it. I used to be able to do it; don't know if something changed or if I'm doing it the wrong way! Something really basic but it is puzzling me and I couldn't find any way that worked in here.


texts:
I have trajectory data in LineStringM format and visualize relevant results via Leaflet on street map tiles. I was wondering whether it is possible to visualize spatiotemporal plots on OSM tiles.
I understand that there are commercial software solutions providing what I am interested in but I was wondering whether it is possible to achieve similar results in a simple manner without drastically changing my workflow using open-source alternatives. It looks like a similar question has been asked here before though the answers might be outdated.


texts: Adding numeric field to end of text string in annotations using ArcMap

texts:
I am using ArcMap 10.8.1.
I carefully moved thousands of annotations for a network that just contain a feature's name. Now for a report map I need to append a numeric score to the end of each annotation.
Current: Parking Lot: 01
Updated: Parking Lot: 01 - (95)
Now, I know this is an incredibly easy process to update the TextString field, but no matter what I do I cannot get the display to update. I just joined a table to the annotation based on the text string and then used the field calculator to append the numeric value to the end of the text string.
So what I have is an annotation class where the attribute table is correct in The TextString column which shows the updated name, but the display still shows the old name. I've restarted ArcMap and even tried opening it on a coworker's copy of ArcGIS Pro.
Annotations are stored in a GeoDatabase


texts: Clipping aerial image pairs by their overlapping area

texts:
I have a set of aerial images which I am trying to clip by their overlapping areas, for use in a more involved program I am writing.

Is there a command line based tool out there with an easy method for extracting the overlapping areas from a pair of images?  The images themselves do not have georeferencing, but they DO come with a shapefile containing bounding boxes of the image collects, as well as IMU data from image collection.

I know this may be difficult since the images will have slight perspective differences, but I was hoping it would still be possible--especially since I am trying to maintain and utilize these perspective differences for something else.


texts:
I am using ArcGIS 10.3.1 and editing a file GDB and a string field.

I have several fields that are similar to the following, 34-01, 23-14, 87-09, 23-12...

I am trying to select for fields that have a "0" as the fourth character.

Is there a easy way to do this?


texts: The number of digit of decimal point when you calculate area of polygon in QGIS1.8.0

texts:
I would like to change the number of digit of decimal point when I calculate area.

For example
Now: 32.456677
Ideal: 32.46

I try to adjust "output field width" and "Precision"
These don't work well.

Let me know how to adjust the number of digit of decimal points.


texts: Get the geometry from the lasso selection tool

texts:
I'm building an Add-in solution in ArcGIS 10 that contains a custom tool. I'd like to borrow the functionality of the lasso selection tool without reinventing the wheel. Is there a way to intercept the geometry of the lasso selection area from the mouse up event? I want to pass the geometry of the lasso area to a spatial filter.

I know how to get control of the lasso tool with the following code in my custom tool class:

private ITool m_SelectTool;
protected override void OnActivate()
        {
            UID uid = new UID();
            uid.Value = "esriArcMapUI.SelectByLassoTool";
            ICommandItem commandItem = ArcMap.Application.Document.CommandBars.Find(uid, false, false);
            ICommand command = commandItem.Command;
            m_SelectTool  = (ITool)command;
        }


texts: How do I add "padding" to features when applying d.out.file to map display on GRASS?

texts:
Thing is this. I've done up a script to automate some spatial interpolation and eventually output the interpolated raster and contour lines as an image of the map display. I'm doing this using d.out.file.

The issue now is, the features (mainly an island + the contours) are rendered right to the edge of the output image. I've tried using d.zoom to zoom out prior to the d.out.file command but to no avail.

Does any GRASS expert here know how I can place a "padding" so that my features don't end up right on the edge of the image? Post-processing won't cut it either as there are grids too. Also, it has to be done programmatically and not via the GUI.

edit: one more thing, I've got a raster mask in place so I'm not sure if that is causing the issue.

Thanks in advance!


texts: Unable to show geometry in virtual layer qgis

texts:
I'm trying to create a virtual layer starting from a municipality layer and from a  table of attributes, both stored in a Postgresql DB.

So I've done a join like:

*SELECT munic, country_of_origin, SUM (presence), munic.NOME
FROM ps JOIN munic ON munic.pro_com = ps.pro_com
WHERE country_of_origin = 'Germany'
GROUP BY munic.NOME
ORDER BY munic.NOME DESC*

I've set the geometry fields in the mask, I see the join in the attribute table, but I can't see the map

Any ideas?


texts: How to show a moving progress bar in QGIS

texts:
I want to sum values and during do it, a window shows progress percentage.
I write below code but when I run it in qgis python console, qgis crashes. what is error in this code?

from qgis.PyQt.QtGui import *
from qgis.PyQt.QtWidgets import *

def progdialog(progress):
    dialog = QProgressDialog()
    dialog.setWindowTitle("Progress")
    dialog.setLabelText("text")

    bar = QProgressBar(dialog)
    bar.setTextVisible(True)
    bar.setValue(progress)
    dialog.setBar(bar)
    dialog.setMinimumWidth(300)
    dialog.show()

def calc(x, y):
    sum = 0
    progress = 0
    for i in range(x):
        for j in range(y):
            k = i + j
            sum =+ k
            progress =+1
            progdialog(progress)
    print sum

calc(10000, 2000)


texts: How to get the exact coordinates of Wolfcamp aquifer locations?

texts: I currently working on the data set "Wolfcamp Aquifer", but the coordinates in the data set is just a relative coordinates. So, is there anyone who has the original coordinates?


texts: Looking for opensource projects for mobile offline trip planning for public transit and walking routes

texts:
I'm looking for initiatives pursuing mobile offline trip planning using a combination of public transit data and walking-routes. Thus far, without much success.

Osmand, seems to be closest, but the part for offline tripplanning isn't completely funded yet. I may contribute, but haven't got a clue about its viability.

I'm looking for any projects pursuing this, or papers (comfortable with reading the academic math-stuff) outlining this. Also roadblocks preventing this (or rather: making quick implementation unlikely), etc. welcome of course.

I'm 'only' looking to do offline trip planning within one city at a time, which obviously brings down complexity lot.

P.S: adapted question from question posted on opentripplanner-forum:
https://groups.google.com/forum/?fromgroups=#!topic/opentripplanner-users/m6PGVRBA0eU


texts: What is the difference between Ground Sampling Distance (GSD) and spatial resolution?

texts:
In research papers on remote sensing, I see some authors referring to spatial resolution (expressed as pixel/distance) and others to Ground Sampling Distance (GSD) (also expressed as pixel/distance). Why would you use one over the other, and what is the actual difference between GSD and spatial resolution?
I found this explanation by Pix4D, but as far as I can tell, GSD and spatial resolution can be used interchangeably.


texts: How to integrate SAGA Python Modules into ArcGIS Desktop?

texts:
I have been having a look at SAGA (http://www.saga-gis.org/en/index.html) and it seems to have some nice tools for geoprocessing.

I was wondering if anyone out there has had any luck bringing SAGA Modules (python) into ArcGIS Desktop and ArcToolbox.

If so what steps did you use?


texts:
I’m looking for the best way to merge text features into one and wrap the text.

Here is what I would like to do in as few steps as possible:

• Select features to be joined
• Merge features attributes (text) (I've been using concatenate)
• Replace "," with "]" since concatenate adds commas ("]" is set to wrap text)
• Set text order ex: 175]+7.5].75”
• Delete duplicate features

There are thousands of features for this to be done so it gets a little tedious to merge the features, add brackets for wrapping and then delete the duplicates. Would there be a way to do this all in one step? Would this be possible to do with the Processing Modeler or Python? (I'm not very fluent with either yet)

Here is an example of how I've been merging features.

Select features to be merged

Concatenate features + replace commas with brackets + set proper order

Then delete the duplicates


texts: Uploading MVT to MapBox result in error "Tile exceeds maximum size of 500k at z5..."

texts:
I have just recently started to use MapBox. I am trying to convert a huge GeoJSON file (1.28GB) into MVT for uploading into a MapBox. I am using ogr2ogr and this command-line:

ogr2ogr -progress --config GDAL_CACHEMAX 2048 -f MVT -s_srs EPSG:4326 -t_srs EPSG:3857 grid.mbtiles grid.geojson

This creates grid.mbtiles and after uploading into Mapbox, it reports an error:

  Upload error:

  Tile exceeds maximum size of 500k at z5. Reduce the detail of data at
  this zoom level or omit it by adjusting your minzoom

I have tried to convert using the next parameters

-dsco SIMPLIFICATION=5 -dsco MINZOOM=5 -dsco MAXZOOM=10

but I still get the very same error.

Could you help me how to solve this issue or give me some advice on how to properly create a vector tileset? How could I determine that the size of particular tiles exceeds this limit before uploading it into a Mapbox studio?


texts: SQL Server STpointfromText comes from PHP into javascript as garbage

texts:
Sorry for the sloppy title.
Here's the scenario - I have a table in SQL Server 2008 R2
It has a geographyColumn, created thusly:

add geographyColumn as geography::STPointFromText('POINT('+convert(varchar(20),x_coord)+'
'+convert(varchar(20),y_coord)+')',4326)

When I viewed it in the table, it looks fine:

POINT (-116.477654 33.764695)

I have some javascript, and a PHP script, to grab points from said table. I can get any text field without a problem. But I want to get the X,Y vales from the geography column.  Instead of a nicely formed "X,Y" or even a long string of numbers, I get this mess:

����zU���f�C@

My php is the following:

<?php $q=$_GET["q"];
$serverName = "ServerIP";
$connectionInfo = array( "Database"=>"geo_spat2", "UID"=>"u", "PWD"=>"p");
$con = sqlsrv_connect($serverName, $connectionInfo);
if (!$con){
   die( print_r( sqlsrv_errors(), true));
};
$sql="SELECT * FROM EMPTY_REPS WHERE Division_label = '".$q."'";
$result = sqlsrv_query($con, $sql);
while($row = sqlsrv_fetch_array($result)) {
  echo $row['geographyColumn'] . ", ";
}
sqlsrv_close($con);
?>

Am I missing something simple?
Thanks,
Evan


texts: Are there any good libraries for accessing GRIB2 data in .Net?  Specifically, I want to be able to quickly pull subsets of weather data from a grib2 file based on time and location.


texts: Display palette (coloramp) imagery in GeoServer

texts:
I have a data set of palette imagery. These are basically RGB imagery that have been converted to 1 band color map imagery using ArcGIS.

I am loading the images to GeoServer as GeoTIFFs but nothing is being displayed in OpenLayers, it is just blank white page. I was wondering if this has to do with their style file, or am I missing something else?


texts: Using Geoserver 2.1.3 I do delete sld using rest api with manually delete the sld file. But when I issue a POST new sld with the same name I got an error "SLD file test.sld already exists." Event there is no test.sld in the styles folder.


texts: How do you update files from OSMDownloader?

texts: After adding the .osm file to your project using the OSMDownloader plugin can the labels or shapes be updated? In other words, when an edit is made in OpenStreetMap and the edit happens to be in your area of interest (i.e. the area you downloaded with the plugin), how do you get the changes to be reflected in your project without having to re-download?


texts:
I can't understand what is wrong in my syntax, or if this is specifically an issue with QGIS but I am hoping someone can explain why I am still receiving a syntax error for this.

I am trying to delete rows based on certain words found in social media posts as they do not apply to my research focus, namely dog-walking.

Here is my SQL query;

DELETE FROM dog_posts_insta
WHERE ((text LIKE '%brew%') OR (text LIKE '%beer%'));

For which I continuously receive a syntax error for near "DELETE"

What am I missing?

My version of QGIS: 3.4 Madeira
Here is a sample of my table as requested:


texts: I want to label a layer with 2 columns/attributes. It´s not a problem in Mapinfo, but can it be done in QGIS as well?


texts: Creating an intersection (mix of offset and buffer) of 5m on each polygon feature in QGIS

texts:
I'm trying to reproduce the steps needed to aim this goal: Given a layer of polygons, I'm want to retrieve for each polygon which is at a distance less than 50m to a 2nd polygonal layer the only 5m band in the polygon facing the items of the second layer.
A drawing may be a better way to explain myself. I want to build the red parts (the 5m band of the orange layer polygons facing the green ones, when they're at a distance less than 50m).
How would you do that with QGIS?


texts: How do I find whether a line is intersecting a circles with small radius on earth?

texts:
Algorithm mentioned in how-to-pick-up-line-segments-which-is-inside-or-intersect-a-circle seems to be working fine for circle with radius > 50 miles.

But when coming to circles with radius < 50 miles it's not working properly. Please let me know an efficient mechanism to handle all the cases.


texts: What is Percentage of Spatial Connectivity - Spatial Weight Matrix

texts:
Generate Spatial Weight Matrix tool in ArcMap and ArcGIS Pro has an output item as "Percentage of Spatial Connectivity". I am wondering how this is calculated and what it represents.
This is an example run output:
Spatial Weights Matrix Summary
Number of Features:                  34
Percentage of Spatial Connectivity:  14.71
Average Number of Neighbors:         5.00
Minimum Number of Neighbors:         2
Maximum Number of Neighbors:         9


texts: I have a problem with update my ArcGIS license manager from ver. 10.5 to 10.6. LM runs on separate virtual machine. My Desktop cannot connect to LM server after update. I've tried to update LM on 2012R2 Standard (x64).


texts: In osm2pgsql how is the planet_osm_roads table populated?

texts:
I'm using osm2pgsql to populate a North America OSM database from planet.osm.

I like how simple the resulting PostGIS schema is, but I'm a bit confused by its contents.

I have the two tables planet_osm_roads and planet_osm_line. I assumed that planet_osm_roads contained all the highway data, but then discovered that it only contains a fraction of the roads in some of the 'highway' types.

For example, planet_osm_line has over 12,000,000 lines with a 'highway' value of 'residential'. However, planet_osm_roads has only 2,000. By comparison both tables have the same number of 'motorway' features. Looking at the data I can't tell why some features would be in the one table but not the other.

I'd really like to see the rules the osm2pgsql uses to direct line features into each table.

I tried looking through the osm2pgsql source code to shed light on the mystery but so far I've not found what I'm after.

Any suggestions very welcome.


texts: Blending DEM hillshade with QMS basemaps

texts:
I would to keep the features of the ESRI Gray (Light) basemap crisp and clear on top of a hillshaded DEM in QGIS.  So far, my attempts have not been very successful.
The ESRI Gray basemap shows main roads in white, and darkened built-up areas on a light gray background.  This is the aesthetic I would like to preserve:

However, overlaying the basemap by blending with "MULTIPLY" on top of a hillshaded DEM tif with a z-factor of 0.00001000 (so as to not create too much darkness and contrast in the hillshade) still really obscures the basemap on top:

I tried making the DEM more transparent, but that fades the landscape features too much, which I would also like to preserve.
Is there a way to get rid of the "greyness" of the DEM that overshadows the features of the basemap, or somehow to keep the sharp features of the basemap to contrast with the DEM?


texts: Organizing PostGIS data of a large database for fast dynamic user searching

texts:
I have a large amount of GIS data, roughly 30 layers, each with 50,000 to 100,000 (Geo)records. Some of this data is GPS collected points, others are POI, some is government provided data (roads, rivers, water sheds, etc),

So far I'll either be using mapfish, geoserver, or feautureserver with jquery/ajax backend.

How should I organize my GIS Data to allow for fast and dynamic webside client searching? as well as feel free to provide your experience from projects especially working on large data sets with heavy user quering


texts: Iteratively unioning pairs of shapefiles

texts:
A simple (I think) Model Builder question. I have 25 field sites, for each of which I have mapped land use as separate layers for different years. In one folder I have veg layers for 1994, named tara{site number}veg94 (eg., tara01veg94.shp, tara02veg94.shp, through to tara25veg94.shp).  In the other folder I have the same for 2007, similarly named by site number (tara01veg07.shp, tara02veg07.shp ... tara25veg07.shp).

I am wanting to union the layers for each site for successive years together using a model, for example tara01veg94.shp unioned with tara01veg07.shp to create tara01veg9407 and so on. I presume that I need to use an iterator to specify the site number for each run of the union function and then insert the %value% created by the iterator into the input and output file names. The problem I'm having is that I need to iterate the same site number for two input files simultaneously, so I figured the "for" iterator was the one to use, but the union function doesn't let me put %Value% into the input fields (though it accepts it in the output field ok).  Any help would be greatly appreciated... thanks.


texts:
I've an array of coordinates that represent a polygon, inspired from this official example I've modified "Fly to Bern" for my aims.

<button id="testID">Test</button>
<script type="text/javascript">

var polygon = [
    [
        13.147888181763385,
        38.696496657431524
    ],
    [
        13.147888181763385,
        38.711500431313084
    ],
    [
        13.16196441467133,
        38.72248334089265
    ],
    [
        13.185997007441344,
        38.7267688965631
    ],
    [
        13.204193113396066,
        38.71069673752169
    ],
    [
        13.197326658318405,
        38.698104355211626
    ],
    [
        13.183593748163982,
        38.68470576954256
    ],
    [
        13.153381345825155,
        38.68470576954256
    ],
    [
        13.151321409301318,
        38.688993589981564
    ],
    [
        13.147888181763385,
        38.696496657431524
    ]
];
var boundaryBox = new ol.extent.boundingExtent(polygon);
console.log(boundaryBox);

function onClick(id, callback) {
  document.getElementById(id).addEventListener('click', callback);
};

function flyTo(location, done) {
  var duration = 2000;
  var zoom = 5;
  var parts = 2;
  var called = false;
  function callback(complete) {
    --parts;
    if (called) {
      return;
    }
    if (parts === 0 || !complete) {
      called = true;
      done(complete);
    }
  }
  view.animate({
    center: boundaryBox,
    duration: duration
  }, callback);
  view.animate({
    zoom: zoom - 1,
    duration: duration / 2
  }, {
    zoom: zoom,
    duration: duration / 2
  }, callback);
};

onClick('testID', function() {
  flyTo(boundaryBox, function() {});
});

  var map = new ol.Map({
    target: 'map',
    layers: [
      new ol.layer.Tile({
        source: new ol.source.OSM()
      })
    ],
    view: new ol.View({
      center: ol.proj.fromLonLat([37.41, 8.82]),
      zoom: 4
    })
  });
</script>

When I click on the button happen that the zoom is on coordinates "0.0, 0.0" instead of the boundary box's coordinates and there aren't errors in console. Probably I've misunderstand something, but what?


texts: How to add documentation for QGIS Processing algorithms?

texts:
In the Processing modeler, many of the algorithm help tabs say:

  Sorry, no help is available for this algorithm.

In the QGIS sourcecode, the processing algorithms seem to reside in:

python/plugins/processing/algs

What general pattern should I follow to add documentation for algorithms (e.g. file naming convention, location, etc)?


texts: Add a popup to a Vector LineString in OpenLayers?

texts:
Is it possible to add a popup (without a marker icon) to a Vector (LineString)?

This is my code, if possible, how to change it to add the popup? I need the popup to display some additional information to the road segement.

<script>
var ytrajectory=[];
ytrajectory.push('LINESTRING (11.1119339 46.0772706, 11.1117177 46.077171, 11.111613 46.0770837, 11.1115376 46.07702)');
ytrajectory.push('LINESTRING (11.1115376 46.07702, 11.1114585 46.0768611, 11.1114219 46.0767673)');

...
</script>
function showReconstructedTrajectoryy() {
var vectorLayer = new OpenLayers.Layer.Vector(
"Reconstructed trajectoryy", {
style : reconstructedTrajectoryStyleYellow
});

var parser = new OpenLayers.Format.WKT();
var src = new OpenLayers.Projection('EPSG:4326');
var dest = new OpenLayers.Projection('EPSG:900913');

var s, wktLine, features, reprojectedFeatures;
for (s = 0; s < ytrajectory.length; s++) {
wktLine = ytrajectory[s];
features = parser.read(wktLine);
reprojectedFeatures = features.geometry.transform(src, dest);
vectorLayer.addFeatures(new OpenLayers.Feature.Vector(
reprojectedFeatures));
}

osmmap.addLayer(vectorLayer);
osmmap.zoomToExtent(vectorLayer.getDataExtent());
}


texts: How to calculate the bounding box in projected coordinates?

texts:
I'm using GDAL in Python to reproject rasters.  I'd like to be able to take an existing raster, and project it into a new raster to my projection and output pixel size of choice.

I'm worried about how to calculate the bounding box in projected coordinates that I will in turn base the geotransform off of.  gdal has a function called gdal.AutoCreateWarpedVRT that apparently creates a memory raster that I can then directly read the geotransform from.  But I sometimes have VERY large rasters that don't even fit in main memory, I'd hate to create a memory based raster just to get the geotransform.

Is it sufficient to transform the vertices of the original bounding box into the projected coordinate system and then somehow make that box a square?  (I'm sure the box will be distorted.)


texts: Calculate proportions of point counts per polygon based on threshold values

texts:
I have an attribute table containing "Joint_Count" and Raster ID (A1, A2, etc), where Join_Count refers to the sum of points per raster ID.

If the sum of all Join_Count represents 100%, how do I retrieve the threshold values to represent 20% and 30% of all values?

FID   Shape    Join_Count        Y_Index     X_Index
0     Polygon      0                   A           1
1     Polygon      1                   A           2
2     Polygon      2                   A           3
3     Polygon      3                   A           4
4     Polygon      4                   A           5
5     Polygon      5                   A           5
6     Polygon      6                   A           5
7     Polygon      7                   A           5
8     Polygon      8                   A           5
9     Polygon      9                   A           5
10    Polygon      10                  A           5

The sum of all Join_Count is 55, corresponding to 100%.

20% and 30% of all valus are:

55*0.2 = 11 and 55*0.3 = 16.5.

Next all Join_Count cells with a sum < 11 (<20%) should be coloured "green", sum values between 11 and 16 (20 - 30%) "yellow", > 16 (>30%) red.

In the above example, this corresponds to Join_Count values < 4 (<20%), between 4 and 5 (20 - 30%), and >5 (>30%).

Any suggestions how this can be done? I'm using ArcGIS 10.3.

I tried to add a new field and run cumsum(!Joint_Count!) in a first step, but I get an unknown error. Also searched through StackExchange, but couldn't find the solution.


texts: Time taken to build the routing network in the pgrouting

texts:
I am using:

PostgreSQL 9.6
pgRouting  2.6.2
PostGIS    2.5

I have line geometry shapefile of count 700K.
When I build the network using  pgr_createTopology in the new DB, network created within 1 hrs.
But When I tried to build the network in actual DB which is a size of 20GB, it tooks to days build.
Before I ran the query pgr_createTopology  in actual DB, I ran vacuum full & checked any temp file present for db in statistics tab.
This is actual settings taken from pgtune:
max_connections = 300
shared_buffers = 512MB
effective_cache_size = 96GB
maintenance_work_mem = 2047MB
checkpoint_completion_target = 0.9
wal_buffers = 16MB
default_statistics_target = 100
random_page_cost = 4
work_mem = 74274kB
min_wal_size = 2GB
max_wal_size = 8GB

I increased the work_mem to 742740kB
Still it took 12 hrs for building 100K records


texts: I try to load views from PostGIS in QGIS 2.18 and I have discovered that QGIS does not recognize the id field of the view and must be set manually. How could I correct this


texts: What happens when you press cancel on a QgsTask?

texts:
I'm trying to get the "cancel" to work for my QgsTask, I have created a simple example to illustrate, When I press the "x" i want it to print reslut: None values: False, hence it just finishes without printing anything. If I let it finish (wait for 20 seconds) its prints reslut: None values: True

import time
def test_task(task):
 for i in range(20):
  if task.isCanceled():
   return False
  else:
   task.setProgress(i*5)
   time.sleep(1)
 return True

def end_method(reslut, values):
 print('reslut: {r} values: {v}'.format(r=reslut, v=values))

task = QgsTask.fromFunction('Test task', test_task, on_finished=end_method)

QgsApplication.taskManager().addTask(task)

Any thoughts on what I shall use instead of task.isCanceled()? (I'm using QGIS3.2)


texts: ArcGIS Geoprocesing Service using map service data for input

texts:
I'd like to create a geoprocessing service in ArcGIS Server (10.2) which uses an existing map service as input.  I have been trying to do this with a ModelBuilder model, but it's copying the existing map service data to my GP service.

More precisely, what I am doing is using a NRCS soil map unit layer and intersecting with a farm field boundary layer.  Both the soil layer and the farm field layer get published with another process.

None of the ESRI examples of GP services show using existing map services as inputs, unless I am missing something.


texts: ArcGIS server Geoprocessing job failed immediately after submit

texts:
I have a few geoprocessing services published onto a production ArcGIS server but after a series of Windows 2012R2 and ArcGIS server updates, the services stopped function.

When I submit a job, these are the messages I would get:

esriJobMessageTypeInformative: Submitted.
esriJobMessageTypeInformative: Executing...
esriJobMessageTypeError:
esriJobMessageTypeError: Failed.

To test and see if there is something wrong with the server, I published an empty python script. When I submit a job to that service, I get the same set of message. If it had function normally, there wouldn't be any "esriJobMessageTypeError". Its as if the code did not execute at all, stopped at the first line.

Does anyone have similar problems after updates and any fixes?

UPDATE
Attempting to built my service from scratch step by step and I think I have located where the problem started.

Starting from the beginning,

Published an empty python script to a test python service and confirmed that it is running normally.
Update the script with some imports. Service is executing normally.
Update the service with one input parameter. The service failed with the above 4 lines of esriJobMessag. However server log did not record anything

UPDATE
Have test the empty script with a single parameter on another ArcGIS server and the server execute the service with no problem.

So what should I do next? Anything before completely re-install the problem ArcGIS server?


texts:
I'm using QGIS 2.16 but this was happening on 2.14 as well and I hoped updating would fix it but no.

I create a new shapefile, doesn't matter whether it's a point, line or polygon, toggle edit function, and then use the edit add feature function to add a feature to the new shapefile. When right-clicking to finish my new feature the feature attributes window pops up but it is blank and I can't edit the attribute which means I can't add the feature. See the image.

I suspect there is a conflict with some other software on my computer or a previous version hasn't been removed completely because I have been having trouble with browsers and videos as well.  I have tried using the tools of CCleaner to tidy up my registry and redundant files.  I tried updating from 2.14 to 2.16 and I have tried removing 2.16 cleaning everything and reinstalling 2.16 but the problem still exists.  Can anybody offer some advice of what else to try?

Note the other edit tools work such as copy, paste, move and node tool and my work around is to copy and paste a feature then reshape with node tool or move it but this is slow.  None of the add... tools work.


texts: Clipping TIF with shapefile containing multiple polygons using Python?

texts:
I'm looking for a way to clip a tif with a shapefile containing multiple polygons like here.

When looking around all the options I find are really long like here and many similar ones. However I can't believe that there is no one line solution for this. Am I missing something?


texts: Rendering and Labeling Shapefile with PyQGIS

texts:
I've been running into a few problems working with the QGIS API for python. I've gotten more comfortable working in the python console in QGIS  but I run into problems when I try to run the code outside of QGIS.

Basically I want to take a shapefile, label it based on the specified attribute name, and render an image. The code works in QGIS, but does not work outside of of QGIS. So where is my problem coming from?

import sys
import qgis
import PyQt4

from qgis.core import *
from qgis.utils import *
from qgis.gui import *

from PyQt4.QtCore import *
from PyQt4.QtGui import *

#initialize QGIS
QgsApplication.setPrefixPath( r"C:\OSGeo4W64\apps\qgis", True )

QgsApplication.initQgis()

#Add layer to instance
file = QgsVectorLayer("Good Shape File", "BMAS", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(file)

#Adjust layer Settings
#Code sample from http://gis.stackexchange.com/questions/77870/how-to-label-vector-features-programmatically
palyr = QgsPalLayerSettings()
palyr.enabled = True
palyr.fieldName = 'Attribute'
palyr.placement= QgsPalLayerSettings.OverPoint
palyr.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,'8','')
palyr.writeToLayer(file)

if file.isValid():
  print "File is valid."

mapRenderer = iface.mapCanvas().mapRenderer()

lst = [file.id()]
mapRenderer.setLayerSet(lst)

mapRenderer.setLayerSet( lst )
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x,y,w,h)
c.addItem(composerMap)
composerLabel = QgsComposerLabel(c)

composerLabel.adjustSizeToText()
c.addItem(composerLabel)
composerLabel.setItemPosition(20,10)
composerLabel.setItemPosition(20,10, 100, 30)

legend = QgsComposerLegend(c)
legend.model().setLayerSet(mapRenderer.layerSet())
c.addItem(legend)

#set image sizing
dpi = c.printResolution()
dpmm = dpi / 25.4
width = int(dpmm * c.paperWidth())
height = int(dpmm * c.paperHeight())
img = QImage(QSize(width, height), QImage.Format_ARGB32)
img.setDotsPerMeterX(dpmm * 1000)
img.setDotsPerMeterY(dpmm * 1000)
img.fill(0)
imagePainter = QPainter(img)
sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
targetArea = QRectF(0, 0, width, height)

#renders image
c.render(imagePainter, targetArea, sourceArea)
imagePainter.end()
img.save("E:/QGisTestImages/out.png", "png")

I'm able to do the simple rendering example in the python cookbook, so I think my paths are set up correctly.

"Good Shape File" should be replaced with a good path location if you want to run this. And palyr.fieldName = 'Attribute' should be set to a valid field name for that shapefile.

Edit:
I have gotten rid of iface and inserted code for the extent between the mapRenderer initialization and the lst declaration.

mapRenderer = QgsMapRenderer()

rect = file.extent()
mapRenderer.setExtent(rect)
mapRenderer.setLabelingEngine(QgsPalLabeling())
lst = [file.id()]

Edit: I added

app = QgsApplication([], True)

after

QgsApplication.initQgis()

and the code worked.


texts: How can I find whether a feature class has editor tracking set to database time instead of UTC?

texts:
I'm trying to publish a map service with about thirty layers in it. During the Service Editor Analyze step, it gave me a high-severity error 00129:

  Services cannot include data where Editor Tracking is enabled and set to record dates in database time. Services can include data where Editor Tracking is enabled and set to record dates in Coordinated Universal Time (UTC).

Suggest solutions (from Esri) include removing the layers, disabling Editor Tracking for the datasets, or turning off field visibility for editor tracking date fields. I'd be happy to do that, but would love to know which layer is causing my problem -- the error only says the data frame has layers, not specific to a single layer.

Is there a quick way to find what layers have editor tracking set to database times instead of UTC? (ArcPy scripts are an option, even if I need to run it on all the feature classes in the SDE -- I just don't know what command or function would check this characteristic!)


texts: Using the Expression String Builder to ultilise specific data within an attribute table

texts:
I have an issue relating to the 'Expression String Builder'. I have a series of data over a 8 year period. I need to collect that data for only a 5 year period. Once collected I need to find the fishing intensity associated with pipelines. From figure 1 I need Dem11 to Dem15 (The figure shows a few pieces of data).

However I don't need Dem07 to Dem10. With this data I open up properties and create rule-basedsymbology. The increments I require are

less than  0
0 - 5
5 - 20
20 - 50
50 - 200
200 - 500
greater than 500

Being new to QGIS, I thought the logical solution would be seen in figure 2.

This seems to work, but not very accurately.

Additionally, I need to do the same with 2 data sets e.g Dem11 to Dem15 & Nep11 to Nep15 together using the same increments. It doesn't seem to give accurate numbers whatsoever. e.g following the formula seen below Dem 11 to 15 gives me (26103 data points between 0 - 5). When I do the exact same with Nep11 to 15. I get (36051 between 0 - 5). However when I use the same formula as seen in figure 3 to add both together to find the fishing intensity I get (23448) which I not sure how is possible as this is 2655 less than the Dem11 to Dem15 data alone. Can anyone help me with my issue?


texts:
I have two shape layers, one containing 93 points and the other layer containing 15075 points. I am trying to randomly assign a new XY position for each of the 93 points to any of the 15075 point locations in the second layer.

I have tried using the create random points tool, but this only lets me generate random points from the 15075 know points. I would like to move each of my 93 points from the first layer to a random location on the second layer so that I can measure the distance moved.

Does anyone have any ideas?


texts: Removing a duplicated coordinate in multipoints using GeoPandas


texts:
Would anyone please indicate to me why this is not working. I am trying to convert all instances of '02ND' to '2ND' for example where some of the input fields (f1) are already converted and some are not.

In codeblock:

 def calc(f1):
 list=['01ST', '02ND', '03RD', '04TH', '05TH', '06TH', '07TH', '08TH', '09TH']

    for i in list:

        if i in f1:
            f1 = f1.replace('0','')
            return f1.replace('0','')
        else:

            return f1

I simply call the function calc to populate a new field.


texts: Creating animations with goTo and extents using ArcGIS API for JavaScript?

texts:
I am using the ArcGIS API for JavaScript.

Right now, I am trying to change the extent of a view dynamically.

<script>
    var view;
    var new_ext;

    require([
        "esri/views/MapView",
        "esri/WebMap",
        "esri/geometry/Extent",
        "dojo/on",
        "dojo/domReady!"
    ], function (MapView, WebMap, Extent) {
        var webmap = new WebMap({ portalItem: { id: "fedda3...." } });
        view = new MapView({
            map: webmap,
            container: "mapdiv",
            center: [-102.45, 47.75],
            zoom: 10
        });

        // Set the extent on the view
        view.extent = {
            xmin: 17.265,
            ymin: 21.99,
            xmax: 17.314,
            ymax: 21.964,
            spatialReference: { wkid: 4326 }
        };
        new_ext = new Extent();

    });

    function react_to_event() {

        // will be loaded by ajax later on
        new_ext = {
            xmin: 17.065,
            ymin: 21.79,
            xmax: 17.114,
            ymax: 21.764,
            spatialReference: { wkid: 4326 }
        };

        view.goTo(new_ext);
    }
</script>

function react_to_event() is just returning an "undefinied".
What would be the right approach to animate the view to an other extent?


texts: Understanding symbology in IDW interpolation using QGIS 2.18

texts:
Continuously to Understand ramp color interpolation using QGIS 2.12, I used IDW interpolation on those height point without inserting the real values (minimum value is 8 meters, maximum value is 584 meters), and got the default ramp values:

But when I change manually the min\max values to the real ones, I get different interpolation:

I try to understand why the values and the symbology changed?


texts: Intersect automatically updating Shape_Area and Shape_Length in ArcGIS Desktop?

texts:
I am making a model in ModelBuilder in ArcGIS Desktop 10.1.

I am using Intersect tool. I need to calculate the Shape_Area and Shape_Length of my polygons (Projected Coordinate System, EPSG: 5514, S-JTSK_Krovak_East_North).

I noticed that these values are automatically calculated and joined to attribute table when I run the Intersect tool.

Is that true for every Intersect run? Or do I need to add tools: Add field and Calculate field to be sure that my calculations of Area and Length are correct and updated after Intersection?


texts: There seem to plenty of resources out there to help users of CAD software move to using desktop GIS but is the reverse also true.


texts: Making meridians converge to both poles?

texts:
def transformx(lon, lat):
    x = lon*(90-lat)/90
    y = lat
    return x, y

What is the formula to make the meridians converge to both poles like the picture attached?


texts: Calculating total line lengths within polygon?

texts:
I have basically the same conundrum as Calculating line segment lengths within polygon using ArcGIS Desktop? but I'm working in PostGIS, not ArcGIS. Any pointers?

This query seems to like it should work, just to show the sums:

SELECT county, ST_LENGTH(ST_Intersection(counties.the_geom,routes.the_geom))
FROM counties, routes
WHERE ST_Intersects(counties.the_geom, routes.the_geom)

But I'm not confident that I'm getting it right and I can't figure out how to reshape that into an update query.


texts: Reverse Engineering of Vector tiles and Tiles Security

texts:
I am using the following workflow to generate and visualize vector tiles.

Store edited shapefile to PostGIS
Create vector tiles using Tippecanoe from PostGIS data
Serve vector tiles through Tileserver-GL

Now when I am inspecting the resultant map site I found the pbf files in Network section. It contains some hexadecimal values. I am using proprietary dataset. This dataset is confidential. I want to know is this secured method? If there is any way of reverse engineering the PBF files of vector tiles what is that and how can I protect my tiles or data security?


texts: Where's the Plugin called Plugin Reloader

texts: Just installed the QGIS 2.16.1 Nodebo version, while trying to sink my teeth into building plugins , I found that the plugin called "Plugin Reloader" doesn't exist !! Please help


texts: Can I shift a georeferenced map a little more left, right,up or down in QGIS

texts:
After having referenced a map, sometimes the algorithm moves the map some pixels to far in one direction. See picture.
In this example, I would like to move the whole lower image more right so it fits more to the upper.
Is there a easy possibility to do so? I could not find a way.
thanks a lot in advance for help.


texts: How to use a feature set as an output parameter in a Python Toolbox?

texts:
I'm having trouble setting a "Feature Set" output parameter from a .pyt ..
I'm using ArcGIS 10.1 SP1.

Here’s a short script that can be invoked either from a classic toolbox (thru the main() function; one setup with same parameters here), or as a .pyt .. it works as expected when invoked from a classic toolbox, but does not as a .pyt:

import arcpy

class Toolbox(object):
    def __init__(self):

        self.label = “FeatureSet Output Test”
        self.alias = “”
        self.tools = [Test]

class Test(object):

    def __init__(self):

        self.label = “Test1”
        self.description = “Output a few points in a feature set”
        self.canRunInBackground = False

    def getParameterInfo(self):

        params = []

        param0 = arcpy.Parameter(
            displayName = “FS Output”, name = “fsOutput”,
            datatype = “Feature Set”, parameterType = “Derived”, direction = “Output”)
        params.append(param0)

        param1 = arcpy.Parameter(
            displayName = “JSON Output”, name = “strOutput”,
            datatype = “String”, parameterType = “Derived”, direction = “Output”)
        params.append(param1)

        return params

    def execute(self, parameters, messages):

        gdbPath = ‘in_memory’
        featureClassName = ‘points’
        featureClass = ‘{}/{}’.format(gdbPath, featureClassName)

        try:
            arcpy.AddMessage(‘Creating a point feature class with 1 text attribute’)

            if arcpy.Exists(featureClass):
                arcpy.Delete_management(featureClass)

            arcpy.CreateFeatureclass_management(
                out_path=gdbPath, out_name=featureClassName, geometry_type=’POINT’,
                spatial_reference=arcpy.SpatialReference(4326))

            arcpy.AddField_management(
                in_table=featureClass,
                field_name=’Attribute1’, field_type=’TEXT’)

            arcpy.AddMessage(‘Inserting a couple points’)
            with arcpy.da.InsertCursor(featureClass, (‘SHAPE@XY’, ‘Attribute1’)) as cur:
                cur.insertRow(((1, 1), ‘point at 1, 1’))
                cur.insertRow(((2, 2), ‘point at 2, 2’))

            arcpy.AddMessage(‘Setting the output parameters’)
            featureSet = arcpy.FeatureSet(featureClass)

            # Has no effect in either .tbx or .pyt case .. no error, just does nothing.
            parameters[0].value = featureSet

            # Does what's expected for classic .tbx; does nothing
            #   for .pyt .. again, no error, just does nothing.
            #   You can comment this line out for testing the first set method ..
            arcpy.SetParameter(0, featureSet)

            # This works fine in both .tbx and .pyt
            parameters[1].value = featureSet.JSON

        except Exception as e:
            arcpy.AddError(str(e))

def main():
    tool = Test()
    tool.execute(arcpy.GetParameterInfo(), None)

if __name__ == ‘__main__’:
    main()

Not sure what I'm doing wrong .. Any suggestions?


texts: Seeking links to free GIS university courses

texts: Does anyone know of any free Study resources for GIS, cartography or Photogrammetry from recognized universities?


texts: Before I go and pay is there a free source for UK Postcode Boundaries, ideally down to Postcode Sector level ("AA00 0" or even better "AA00 0A")?


texts:
I want to findall tags that match

root.findall("{http://trimble.com/schema/fxl}FeatureDefinitions/{http//trimble.com/schema/fxl}PointFeatureDefinition")

or

"{http://trimble.com/schema/fxl}FeatureDefinitions/{http//trimble.com/schema/fxl}PolygonFeatureDefinition"

or

"{http://trimble.com/schema/fxl}FeatureDefinitions/{http//trimble.com/schema/fxl}PolylineFeatureDefinition"

I need to be able to grab anyone of those three tags, I think the parenthesis is throwing me off. I tried some regex and findall("*") works but after I go further I get no return.


texts: Measured line distance very different in ArcMap than Google Earth

texts: My question is about why Google Earth and ArcGIS have such dramatic differences when measuring a line-string. I have an approximately 30 mile trail in Butte County, Northern California. I measured it in ArcMap using 3D analyst, Add Surface Information and draping it over a DEM and it comes to approximately 28 miles. This is UTM 10N. When I export the line as a KML and bring it into Google Earth and look at the elevation profile, it tells me that the length is 33.3 miles, which is wildly different than expected. Does anybody have any thoughts as to why this is so?


texts: qgis Grouping Analysis - Assign points to subgroups - concavehull not working

texts:
I have a series of polygons I have converted to points.  I want to have qGIS automatically assign them into subgroups of 15 points that are near each other.

So I want something like this to occur:

I found a tool description in archGIS that appears to do exactly what I want, but can't find the equivilant in qgis.  (The tool is here: Grouping Analysis (Spatial Statistics)

I can't figure out which tool in QGIS to use.  I've seen elsewhere to try the concaveHull plugin, but I can't get that to work right. Regardless of what I set the 'Number of neighbors' and 'Find SNN Clusters' values to, it creates just a single shape with only one record.

What am I missing?


texts: How should I use features when iterating a feature class with a cursor?

texts:
I'm iterating through all the features of a feature class (polygons) and I don't think I understand very well what each row is exactly. See example below of what I'm trying to do with arcpy:

DSM_model = os.path.join(working_folder, 'DATA\Data4Optimimty\ESRI_UK_Bluesky_DSM\ESRI_UK_Bluesky_DSM.gdb', 'DSM_25cm_2015')
fc_London_Bldngs_DSM = os.path.join(working_folder, 'DATA\Optimity_Temp.gdb', 'London_Bldngs_DSM')
field_names = ['OBJECTID']

with arcpy.da.SearchCursor(fc_London_Bldngs_DSM, field_names) as cursor:
    for row in cursor:
        DSM_bldng = os.path.join(working_folder, 'DATA\Optimity_Temp.gdb', 'DSM_bldng')
        temp = ExtractByMask (DSM_model, row)
        temp.save(DSM_bldng)

I get this error:

 Traceback (most recent call last):
  File "C:\ESRI\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\Find_Flat_Roofs.py", line 44, in <module>
    temp = ExtractByMask (DSM_model, row)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 7249, in ExtractByMask
    in_mask_data)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\sa\Functions.py", line 7245, in Wrapper
    out_raster)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool
 Failed to execute (FindFlatRoofs).

I believe the second parameter of 'ExtractByMask' is expecting a feature layer, but I don't know what 'row' is.


texts:
I am trying to understand if the Landsat 7 satellite ETM+ orbit has any significant azimuth or if it could be assumed as 0° (as it is not a pointable sensor).
I searched in the web but could not find any good answer so far.
Zenith is supposed to be equal zero as Landsat images are always nadir.

What about the azimuth?

Would knowing the scene center coordinates (they're 46.036323 N and 10.456876 E) be of any help?

The image looks like a sort of isosceles trapezoid, with the major basis towards the North (as I would expect for the Northern Hemisphere). It is also tilted differently from south to North (as I also would expect).

One question arose while I was trying to solve my task. The metadata file of Landsat provides coordinates of the four image corners. Anyway they refer to the black background.

The azimuth is intended to be calculated through these corners or the corners of the "real"-non-background image?

And why do I have to use latitude rather than longitude?

Is it not better to use longitude?

I thought to use (for instance) the right corners, subctracting the longitudes of them (e.g. upper right long-lower right long) and obtain what I was looking for in terms of degrees, does it make sense for you?

Here is an image showing a part of Landsat paths in a Worldwide Reference System (descending orbit) and my Landsat scene. RS is geographic, WGS84. According to this what I should looking at is the non-background image.


texts:
I want to import several spatial csv files (points) into QGIS.
The path to the folder, where the csv files are located is stored in the variable "Folder_Path".
The file names are stored in a list.

It seems to me, that the only way to create a memory layer in this case is with a uri. But the uri in the syntax of QgsVectorLayer(uri, 'Test', 'delimitedtext') has instead of a forward slash "/" a backward slash "\".
Replacing the forward slash with a backward slash doesn't work with my code:

Folder_Path = QFileDialog.getExistingDirectory(caption="Pick Folder 'Import_QGIS'.")

#Folder_Path for uri
Folder_Path_uri = Folder_Path.replace("/", "\")

Result is SyntaxError: EOL while scanning string literal

How can I let Python know, that the "\" is just a character and not an escape?
Or is there a better way to go about adding csv files?


texts: Using "texture" image from Drone2Map for ArcGIS?

texts:
ESRI's  Drone2Map beta produces various outputs; LAS, dem, hillshade, orthomosiac, etc.  One of these is "Texture.jpg" file.  This image is not projected but it looks like as if I should be able to apply to the sides of structures and vertical "faces" of features.

Does anyone know how this can be done in ArcScene and if not in ArcScene, what sofware would allow me to do this and how?

capture 1: UAV data processed in ESRI's Drone2Map and viewed in ArcScene (notice vertical faces are stretched.)

capture 2: a top-down view of the "texture" image created by the Drone2Map process


texts:
I try to measure the length (in meters) of the polyline layer that pass inside the olive polygon layer only, and if it possible export the wanted polyline into new layer:

the attribute table of the polyline is:

All the layers as the same CRS.


texts: Need links on DWG to GDB conversion using FME Workbench

texts:
I am using FME workbench (trial version) to convert .DWG to .GDB (ESRI Geodatabase [file based]). I have mapped the Writer attributes with the exact Writer attributes. In fact, in the test database that we are using, there are only 5 fields with char(255) datatype, so no problems while mapping them.

Now when I run translation, it says 'Translation was SUCCESSFUL'. But also gives following warning:

During translation, some features were read that did not match a reader feature type in the workspace. This can happen if the reader dataset is changed, or a reader feature type removed or renamed.

If I test in ArcCatalog, the destination FeatureSet is created will all shapes, but WITHOUT fields mapped from Reader to Writer. If I preview table, it displays the table with all the rows with only Object_ID, Shape and Shape_Length fields, but not the attributes that were mapped from Reader to Writer.

I am aware that it might not be possible for you give any comments without having a look at the workspace. But can someone provide me some links to tutorials showing .DWG to .GDB conversion using FME workbench.

Thanks in advance


texts: Is there an R equivalent of NEAR in ESRI ArcGIS?

texts:
I have 54 points in one data set A
I have 300+ data points in a dataset B.

I wish to find the distance for each point in dataset-set B, to points in dataset A.

Nearest Neighbor packages appear to only find neighbors within a given data-set, not a second data set.

rdist was close, but I didn't want a matrix, just nearest.--"Given two sets of locations computes the full Euclidean distance matrix among all pairings or a sparse version for points within a fixed threshhold distance."

This looked good, saying "Read two geocoded point sets from Comma Separated Value (CSV) files into R data objects. Assign to each member of first point set the geographically closest point from second set" but is so dated (2010), so something more recent must exist...

Same issue with this

Other packages find distances between points and lines, points and polygons, etc.
viz: http://cran.r-project.org/web/packages/geosphere/geosphere.pdf

The omission of P2P suggests I am missing something simple...


texts: Identifying missing tiles on basemap on ArcGIS Server?

texts:
I have a basemap service on ArcGIS Server 10.1 where after doing a cache update (exploded cache), some of the tiles are missing.  I have added a "missing tiles.png" so that it shows where the missing tiles are.

However, i have a very large area cached at 3 scales.

Using ArcGIS Server or ArcGIS Desktop (Advanced level license) is there a way to find out which tiles are missing and at what scales?

The cache has been created using my own tiling scheme file (created through ArcGIS Desktop).

Ideally I would like to know which tiles from my tiling scheme are missing from my basemap service cache.  Finding them visually will be extremely time consuming.

I've looked in the "Status.gdb" for the cached service, but there doesn't seem to be anything useful in there to tell me where the data is missing.

I've tried doing a cache update to "Recreate empty tiles", but this doesn't seem to work, and as it is a very large dataset it takes days to process.

I was hoping to be able to extract a polygon/number of polygons for the areas of missing data and send those in as the AOI for a cache update, as this will be a much smaller subset of the data.


texts: How can we split a polygon within another polygon in Qgis?

texts: How can I split the green area as a separate polygon fro the grey one in QGIS?


texts: I am joining points to polygons in ArcGIS (10.5), but this is taking a ridiculous amount of time (roughly a week!). The points are climatic data covering the UK - approx 450,000 points. The polygon is a fishnet covering the UK, each cell is 500m^2 (approx 3.8m cells). I don't know how to find the important information that may help to understand why its taking so long, e.g. data formats, volumes etc..


texts: QGIS 3.2 Raster Histogram frequency axis always locked at 1

texts: Question is pretty much in the title: Whenever I go to layer properties and generate a histogram for a raster the frequency axis of the histogram is capped at 1. So for say a multi-band image I see the pixel values for the multiple bands on the x-axis of the histogram and the frequency of those values on the y-axis. My problem is that the range of the y-axis displayed is from 0-1 which is useless. How do I fix this?


texts: How to convert Blob data to string in Arcpy?

texts:
I am using Arcpy to load a table from Excel. One of the fields has a lot of text (up to 2000 characters) and when using arcpy.MakeTableView_management the field is given the type 'BLOB'. Is there any way to either prevent this from happening, i.e. making sure the field type is 'Text', or is there a way of reading the text in python so the data can be used?

I have read BLOB TYPE in Arcmap but it does not answer the question. I am specifically trying to read the data in Python, while the previous question referred to reading it in ArcMap.


texts:
I have installed QGIS (2.14.3-Essen) along with SAGA GIS (Version: 2.2.7) in Ubuntu (14.04 trusty). I am trying to fill the sinks in a DEM in QGIS processing using SAGA provider. I have used both the conventional Fill sinks and Fill sinks xxl (wang & liu) commands. But it gives me the same error.

Following is the content of the log for the command. The error is seen at the last line, where it says that the Filled DEM could not be generated correctly. Could anybody please help me with this error?

Following input parameters were supplied:

DEM: A 30m SRTM DEM
Minimum Slope (Degree): 0.10000 (Also tested with 1.0000 and 0.01000, and that also gave the same error).

Log output:

Algorithm Fill sinks xxl (wang & liu) starting...
ta_preprocessor "Fill Sinks XXL (Wang & Liu)" -ELEV "/tmp/processingd598aea9706848d19dfbd7650d1adaec/dc79aa0c5343429683c3aedae3409e64/clippedbagmati.sgrd" -MINSLOPE 2 -FILLED "/tmp/processingd598aea9706848d19dfbd7650d1adaec/f96c3d2452e642d789f4b534cf63960b/FILLED.sdat"
____________________________

##### ## ##### ##
### ### ## ###
### # ## ## #### # ##
### ##### ## # #####
##### # ## ##### # ##
____________________________

SAGA Version: 2.2.7
____________________________
library path: /usr/lib/x86_64-linux-gnu/saga/
library name: libta_preprocessor
library : Preprocessing
tool : Fill Sinks XXL (Wang & Liu)
author : Copyrights (c) 2007 by Volker Wichmann
processors : 4 [4]
____________________________

Load grid: /tmp/processingd598aea9706848d19dfbd7650d1adaec/dc79aa0c5343429683c3aedae3409e64/clippedbagmati.sgrd...

Parameters

Grid system: 29.444285; 3159x 2645y; 600784.721897x 3000440.861314y
DEM: clippedbagmati
Filled DEM: Filled DEM
Minimum Slope [Degree]: 2.000000

Converting outputs
Loading resulting layers

The following layers were not correctly generated.
Filled DEM
You can check the log messages to find more information about the execution of the algorithm

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


texts: How to get the editor's current layer in ArcPy?

texts:
In ArcObjects, it is possible to retrieve the currently selected editable layer using IEditor3.CurrentTemplate.Layer, or IEditLayers.CurrentLayer.

Is it possible to do the same thing in ArcPy?

In Python add-in extensions, there is an onCurrentLayerChanged method, but no corresponding property as far as I can tell to determine the current layer, which begs another question -- what is that method intended to be used for?


texts: How to make the white background of an RGB raster transparent?

texts:
I would like to show on my map only the reddish pixels of my raster layer.

Said in another way: I would like to transform my white background into transparent background so I can see my other features.

I have the feeling this is possible using the custom transparence tab in the layer property.


texts: For looping folder to batch clip rasters by polygon using python and QGIS?

texts:
I am using python and QGIS 2.0.
I am trying to clip rasters in a folder by one polygon feature.
It's the first time for me using (let's say) "PyQGIS", I was used to arcpy before.
Anyway, I don't get my simple script to work, any suggestion would be very appreciated!

import qgis.core, qgis,utils
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()

CLIP= "C:/Users/unim/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/study_area_foscagno.shp"
INPUT_FOLDER="C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/LE71930282000259EDC00"
OUTPUT= "C:/Users/unim/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/foscagno_pyqgis/"

for RASTER in INPUT_FOLDER.tif
do
    echo "Processing $RASTER"
    gdalwarp -q -cutline CLIP -crop_to_cutline -of GTiff RASTER OUTPUT+ "clip_"+ RASTER
done

QgsApplication.exitQgis()

Below are the improvements I've done since now, not getting the script to work though, but I think I might be getting closer...

import qgis.core, qgis.utils, os, fnmatch
from osgeo import gdal

CLIP= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/study_area_foscagno.shp"
INPUT_FOLDER= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/LE71930282000259EDC00/DNs2Reflectance_LE71930282000259EDC00"
OUTPUT= "C:/Users/unimi/Documents/Umberto/Universita/PhD/Guglielmin/Permafrost/Alta_Valtellina/Landsat_ita/Cloud_mask_AltaValtellina/clip_2_foscagno"

def findRasters (path, filter):
    for root, dirs, files in os.walk(path):
        for file in fnmatch.filter(files, filter):
            yield os.path.join (root, file)

for raster in findRasters (INPUT_FOLDER, '*.tif'):
    print (raster)
    outRaster = OUTPUT + '/clip_' + raster
    cmd = 'gdalwarp -dstnodata 0 -q -cutline CLIP -crop_to_cutline %s %s' % (raster, outRaster)
    os.system (cmd)

I think there might be something wrong in the "gdal" command, as the "print" function does its job properely, but no file is written to the output, neither I receive any error. By the way, it's been hard to fond an easy documentation on gdal coding...


texts: Confusing start times vs. time windows in Network Analyst

texts:
This is the snippet of my ArcPy script that deals with route-making:

TechRoute = arcpy.na.MakeRouteLayer(inNetworkDataset, RouteName, "Minutes",
                                    "FIND_BEST_ORDER", "PRESERVE_FIRST", "USE_TIMEWINDOWS")
field_mapping = "TimeWindowStart TWStart1 #;TimeWindowEnd TWEnd1 #"
arcpy.na.AddLocations(RouteName, "Stops", TechLayer, field_mapping, "1000 Meters", "VisitOrder")
outputLayerFile = os.path.join(workingDir, "opt{}".format(RouteName))

Note that I'm not assigning any value for the start_date_time parameter for na.MakeRouteLayer. I decided to leave it out since I want the route to start at 8am, and in the Esri help page, it says

  The default value for this parameter is 8:00 AM.

Unfortunately, it's not starting anywhere CLOSE to 8am. (And also not starting on today's date, but I am not interested in the date so I don't care.)

As best I can figure it's using the time windows, which (for this particular route) all start at 8am and end at 6pm, and working backwards. This isn't what I expected and also seems a bit silly (wouldn't most users want to start at the beginning of a time window?). For another technician's route without any appointment windows, it's starting at midnight instead of 8am. (Again, not what I would expect!)


texts: Deactivate Snapping for Measure Maptool (mActionMeasure) in QGIS temporarily

texts:
I use the QGIS measure-tool to measure distances. I often would like to deactivate the snapping temporarily while using the measure-tool.

Is it possible to deactivate snapping (for example with a shortcut) while measuring and re-activating is when i digitize again?


texts: Intersection of 4D Trajectory with 4D Polyeder (3D + time)

texts:
I have worldwide small polygons (lat/lon WGS84) for different heights (z-coordinate) for several timestamps (m). I want to store them in a PostgreSQL/PostGIS database as 4D objects (2D polygons + height range + single timestamp).

I have several routes via multiple waypoints in the world with different start/end times. I know the time and height at every waypoint.

I want to get all of the polygons that are intersected by the route in four dimensions (lat/lon/z/m). Ideally with a buffer around each of the dimensions.

Is this possible with a single SQL SELECT with PostGIS? Are there other approaches?


texts: Erdas mosaic producing lots of black space?

texts: I am trying to create an imagery mosaic using erdas imagine. The original files were in a 16 bit DIM Pan-Sharpened format, which produced a completely unusable blank raster (likely an issue with color balancing). I used ArcMap to turn every image into a 8 bit JP2, but that produced images with lots of black space in order to create rectangles, rather than rhombi. The problem being, is that erdas did not clip out any of the empty space when it was mosaicking. The only possible solution I can think of would be to produce the seamlines in Arc, then use those same seamlines to clip with in erdas in a shapefile. When I tried to load the shapefile into the erdas mosaic tool, it produced an error. Has anybody else dealt with this problem, and is there an easier way that simply clipping and mosaicking one by one in ArcGIS?


texts: How to draw a line between two places with via-points?

texts:
I have prepared a GPS-tracking dataset, currently consisting of one table with stopping points and one table with the raw GPS-points. The first table contains a column for any stopping point showing the id's of the grouped points.

Now, I want to connect the stopping points with lines, not only direct connections (st_makeline(geom, geom)), but with via points from the raw point table. Therefore, I tried something like

st_makeline(array[stoppingpoint, viapoints_lying_between_these, next_stoppingpoint])

Unfortunately, this query connects all via-points with the stopping point with one line, as you see here:

The query looks as follows:

    st_makeline(stops.the_geom, raw_points.the_geom, ARRAY[lead(stops.the_geom) OVER (ORDER BY stops.group_id)])AS lines,
    stops.group_id, row_number() OVER (ORDER BY stops.group_id) AS line_id
from // semantic: where raw_points.id lies between maximum_id_from_first_stop and minimum_id_from_next_stop

When I manually put the via points and the two stopping points in one array within the st_makeline-function, it works in terms of drawing one line.

Thank you for any advice!

UPDATE:

okay, i did not solve the problem yet. So, here is the table for better understanding: What i need are lines from each stopping point (marked by column group_type), which are centroids of grouped points, to the next stopping point VIA the intermediate points. In this case, one line would start at the point with group_id 38, end at group_id 42 and be connected via points 498, 499 and 500.

UPDATE:

I finally found a solution. First i merged all the consecutive waypoints (id 38-42 in the lower picture) with a window function giving new group_nr's:

   sum(group_flag) OVER (ORDER BY id) AS group_nr
   FROM ( SELECT
          CASE WHEN lag(group_type) OVER (ORDER BY id) = group_type THEN
            ELSE 1
                END AS group_flag from table) t1

Then i make lines from each group of waypoint (multipoints) and add the previous and the next stopping-point-geometry to each line:

    -- 1. previous stop:
st_addpoint(st_linefrommultipoint(groups.the_geom), lag(groups.the_geom) OVER (ORDER BY groups.group_id), 0)

    -- 2. next stop:
st_addpoint(t1.new_geom, lead(t1.the_geom) OVER (ORDER BY t1.group_id)


texts:
I've a layer who contain many parameters like number. I try to "navigate" in this layer for make a specific action.

The schematic structure of my layer look like :

object : [a] - item (1) - [o] - name (Name layer) - [o] - source (14) - [o] - number (14)

for take the name of my layer I'wrote (return "Name layer") :

map.getLayers().item(1).get('name')

but if I want to make an action on the layer "Name layer" and the number "14", how can I do that ?

I tried : map.getLayers().item(1).get('name').item(14).get('number');

But it's don't that.


texts: Method to determine a person's city/state using location coordinates?

texts: I am working on data-mining project where I have to map a bunch of tweets to US states. I can get the lat/long coordinates from the tweet json, but I am not sure if it is possible to map those coordinates to a state? Does anybody have idea/approach to how to approach the problem? Is there any service which I can use for my project? or a database which maps coordinates to US states?


texts: Elegant ways to express coordinates and zoom level in urls

texts:
To make sharing map locations more natural I'm looking at updating a map service to sync the browser url with the current center point coordinate and zoom level the user is viewing.

At present the service is using uglier url parameters e.g.

map?lat=-43.273331&lon=171.960411&zoom=12

What other friendly formats are in use out there already or could you suggest?

Some methods I've seen in the wild include:

http://maps.gsi.go.jp/#16/40.649778/140.268438
https://www.google.com/maps/place/Amsterdam,+Netherlands/@52.3545362,4.7638767,11z/

I realize some of you will suggest using a bounding box instead of a single center point, but in this case different topo imagery is shown at different zoom levels so I'd like to preserve the level of detail the current user is looking at when shared. Topo250's and topo50's in this case.

I feel the url parameter approach is uglier as they are longer and potentially mixed in with other page url parameters.


texts: PyQGIS standalone script error message for gdal_merge: "Incorrect parameter value for INPUT"

texts:
In my standalone PyQGIS script I would like to merge two raster layer. "zero" and "raster". Both variables are definitly QgsRasterLayer objects. This is the function:

def merge(self):
    zero = self.setToZero()
    raster = self.getRasterLayer()
    outDir = os.path.join(QgsProject().instance().homePath(), "Verschmelzung.tiff")
    param = {"input": [zero, raster],
             "output": outDir}
    merge = processing.run('gdal:merge', param)
    rlayer = QgsRasterLayer(merge['output'], 'mergedRaster')
    return rlayer

An that is the error message:

_core.QgsProcessingException: Unable to execute algorithm Incorrect parameter value for INPUT  Traceback (most recent call last):   File "C:/Users/Jan/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\gebeude_betroffenheit\Gebeude_Betroffenheit.py", line 347, in run
    self.zonalStatistics()   File "C:/Users/Jan/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\gebeude_betroffenheit\Gebeude_Betroffenheit.py", line 315, in zonalStatistics
    raster = self.focalStatistics()   File "C:/Users/Jan/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\gebeude_betroffenheit\Gebeude_Betroffenheit.py", line 297, in focalStatistics
    raster = self.merge()   File "C:/Users/Jan/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\gebeude_betroffenheit\Gebeude_Betroffenheit.py", line 292, in merge
    merge = processing.run('gdal:merge', param)   File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python/plugins\processing\tools\general.py", line 105, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)   File "C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python/plugins\processing\core\Processing.py", line 139, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException:
    Unable to execute algorithm Incorrect parameter value for INPUT

According to this link it wants a list of rasters, and that is what I'm giving no.

I don't see the problem…

I'm using QGIS 3.4


texts:
What I try to achieve is quite simple. I want to filter my vector layer by a query in the features.
What I found out is that vector layer doesn't have updateParams() where I could do something like:

'layerDefs': featureName1+operatorName1+checkValue1

So is there an other way to achieve this? Cause now I`m kinda clueless on how to filter a vector layer.
Every way is welcome, the vector is already made so it's not a new vector layer.

I`m working with OpenLayers3 and ArcGIS featureservice.

I`m also not sure how to redraw the vector after it's filtered yet.
All of my vectors are saved in an array.

I tried everything but sadly i couldn't get it to work. this is what i tried in the end. Last try is underneath with help of the suggestion from @Alexis_A

Okay after some more testing i found the problem, for some reason it gives the name of my variable and not the value that contains my variable.
So in other words this bit of code:

var fn = featureName1;
console.log(fn);
var filtered = features.filter(function(feature) {
    console.log(feature);
    console.log(feature.values_.fn.value);
    return feature.values_.featureName1 == 2;
});

will littarly give this as console log:

  categorie_melding

  Object { pendingRemovals_: Object, dispatching_: Object, listeners_:
  Object, revision_: 2, ol_uid: 90, values_: Object, id_: 158816,
  geometryName_: "geometry", style_: null, styleFunction_: undefined, 2
  meer… }

  TypeError: feature.values_.fn is undefined[Meer info]

So at first it knows what the variable should be, but when i`m in the filter function it's just clueless. Anyone any idea on how to fix this?

And why is filtering in OpenLayers3 so hard i just can't understand that they dont have an easier way to do it.


texts: Vector -> Geometry Tools -> Extract Nodes equivalent in PyQGIS API

texts:
My question is fairly straight-forward.

Would anyone know the PyQGIS equivalent of the Extract Nodes tool in QGIS? I don't seem to find it in the API documentation. If you do, kindly share the documentation too.


texts: Converting 2D point cloud to polygon using ArcMap

texts:
There was similar questions, but not exactly what I need, so I am asking here.

Is there a way to create polygons from these points in the picture?

The result I want to achieve in this picture >>

The points in the picture is GPS points from transport working in the field. From those points I want to have polygon geometries for future calculations of area and boundaries.

Zoomed out version of my GPS points in 1:1,500,000


texts:
What is that borders in the water on OpenStreetMap maps?

Is this some of maritime zones?

Is this "territorial waters", i.e. 12 nautical miles away from borders?

What is between islands?

Where/how to obtain this information?


texts: Display mouse position from graph on a map

texts:
i need to develop a web application which shows a map, that displays a GPS Route.
So far I can do this using Openlayers or ArcGIS Server.

Now I want to add an interactive graph below the map. In this graph I want to display for example "Speed" and on mouse-over on the Graph the corresponding location in the map should be displayed.

Does anyone know an existing library or can give me recommendations on which architecture to use? Or do I need to programm this thing from the ground up on my own?

Thanks!


texts: How to ST_Split features in one table by features in another?

texts:
I need to split polygons (layer 'pol') by both closed and non-closed linestrings (layer 'lin').

Unfortunately I don't get proper results running the following query.

CREATE VIEW splitted_pol AS
SELECT
    g.path[1] as gid,
    g.geom::geometry(polygon, SRID) as geom
FROM
    (SELECT
    (ST_Dump(ST_Split(pol.geom, lin.geom))).*
    FROM pol, lin
) as g;

In my example ST_Split should create six polygons (layer 'splitted_pol').

 Does anybody know how to use ST_Split with QGIS/PostGIS?


texts: How to use dojo or esriRequest object to get pdf that is stream in a byte array

texts:
I have an SOE that creates a pdf and streams it back in the response and it works great.

How can I call it using dojo or the esriRequest object to call the SOE and display the resulting pdf in a new window?

I see the esriRequest has a handleAs that can be set to document, but I fail to see how it is implemented.  I have seen examples of using dojo.io.iframe, but I prefer not to use a deprecated library.

var downloadPdfIframeName = "downloadPdfIframe";
var iframe = dojo.io.iframe.create(downloadPdfIframeName);
dojo.io.iframe.setSrc(iframe, url, true);

I would be nice if I could just use straight javascript and do something like window.open(encodeURIComponent(url)), but the url would be too long after I add all the parameters to the querystring and it would have to be a post request.


texts: PostGIS raster - how to populate pixel values by set of points?

texts:
I have the table "events" with column "point" that contains geometries (points). And I need to fill a raster (already created in table "rasters") using this set of points.
For simplicity BandPixelType is 1BB, default pixel value is 0, and I want it to be 1, if pixel contains at least one point from my table. Raster and geometries have the same SRID (4326).
So, I need something like:
UPDATE rasters SET rast = ST_SetValue(rast, 1, ST_SetSRID(ST_Point(lon, lat), 4326), 1)
but executed with a dataset:
SELECT point FROM events WHERE ...condition...


texts: Showing number of features in Contents pane of ArcGIS Pro

texts:
Is it possible to show the number of features per category in the Contents pane of ArcGIS Pro?

It is simple to do in QGIS, but in ArcGIS Pro I only found the possibility to retrieve the number of features per category within the Symbology pane.


texts: Is ESRI's ArcGIS platform of programs available on Ubuntu? I love open source and if I could I would love to try and migrate to this open source OS.


texts: Trilateration with error account or lack of precision

texts:
I am developing a node.js/javascript application where I would have dozens of beacons in a row (positions known), which gives me their distances to specific object (3 or 4 distances per object). The idea is to trilaterate those distances to infer the coordinates of this object.

I have technical background but I am not an expert in geometry nor cartography. I have check a lot of articles about trilateration. Most interesting:

Trilateration algorithm for n amount of points
Trilateration using 3 latitude and longitude points, and 3 distances

With the distances involved (hundreds of meters), instead of 3D problem, I can simplify to 2D problem, but I am not sure which would be best approach to tackle it. I would need a fast algorithm, since I can have around 12 different position calculations (different objects) every second.

I have different ideas:

Using wikipedia's article is the easiest way. I have a pair of javascript libraries that perform this algorithm (cannot put links due to constraints in stackexchange). Problem: does not take into account error in measurements, so many times equation does not give you a solution (approximate value would be enough...)
Using Mathematica: cannot do that. Is there any other C++ or javascript/nodejs library that performs those calculations?
Using postgis: I think I can create a circular stripe (donut) per beacon distance, using ST_BuildArea. Then intersect those areas and calculate the centroid. Do you think this would be a proper approach using geometry - not geography-? Fast enough for 12 positions/second?


texts: Using Java, I want to add geo tags to existing TIFF file and make it into a GeoTIFF file. Otherwise, for buffered image is it possible to add geotags to make a GeoTIFF file.


texts: How do I select polygons that overlap by a specified percentage of overlap?

texts:
I have two polygon layers. One represents parcels, the other is a hazard zone. How do I return only parcels where more than 30% of the parcel area overlaps the hazard zone? For example a parcel with an area of 100,000 square feet with 35,000 square feet inside the hazard zone should be returned as it's 35% overlap. On the other hand, a parcel of 100,000 square feet with 5,000 square feet inside (overlapping) the hazard area should not be returned as it's only 5% overlap.

I tried st_intersects but, although I get the parcels that overlap the hazard are, I need to take into account the percentage of overlap. How do I compare areas of overlap versus parcel area?

select p.id from parcels p, hazard h where st_intersects(p.geom, h.geom)


texts: layer.setQuery is not a function error from CartoDB?

texts:
I am getting this error:

SELECT * FROM parcels_sa WHERE parcel_id ILIKE '%001089800%'
mapscript.js?1437514849:273

Uncaught TypeError: layer.setQuery is not a function(anonymous function)
@ mapscript.js?1437514849:273x.event.dispatch @ jquery.min.js:5x.event.add.y.handle @ jquery.min.js:5
cartodb.js:5 GET http://dchadney.cartodb.com/api/v2/sql?q=SELECT%20ST_XMin(ST_Extent(the_geo…arcels_sa%20WHERE%20parcel_id%20ILIKE%20%27%25001089800%25%27)%20as%20subq

  net::ERR_EMPTY_RESPONSE

function parcelSearch(layer) {
    var sql = new cartodb.SQL({ user: 'username' });
    var $options = $('#searchParcel');
    $options.click(function(e) {
        var streetNam = document.getElementById('parcelid').value;
        if(streetNam === "") {
            alert("Please enter the parcel ID");
        } else {
            query = "SELECT * FROM parcels_sa WHERE parcel_id ILIKE '%" + streetNam + "%'";
            console.log(query);
            layer.getSubLayer(0).show();

            sql.getBounds(query).done(function(bounds) {
                latlon = bounds[0];
               // console.log(latlon);
                //alert(latlon);
                $('.latlongit').val(latlon);
                map.fitBounds(bounds);
            });
        }

        // change the query in the layer to update the map
        layer.setQuery(query);

    });

    var $clear = $('#clear');
    $clear.click(function(e) {
            $('#parcelid').val('');
            query = "SELECT * parcels_sa";
            console.log(query);
            layer.getSubLayer(0).hide();
            sql.getBounds(query).done(function(bounds) {
            latlon = bounds[0];
            map.fitBounds(bounds);
            });
        // change the query in the layer to update the basemap
        layer.setQuery(query);
    });
}


texts: Making one field equal to another in the same Attribute table using Field calculator in QGIS

texts:
I am trying to have a new field with the same content of an existing field.
In ArcGIS I would normally use the Field Calculator to make this instantly.
In QGIS I have followed the same steps (create new field, open field calculator, make it equal to an existing field) but it doesn't work... meaning nothing changes, the new field remains empty...

Any suggestions?


texts: How to Iterate through Bands of Raster using ArcGIS 10.1 ModelBuilder?

texts:
I have a large number of rasters with 3 bands (R,G,B). What I want to do is add all 3 bands of each raster together to end up with a single band image with the summed R,G,B values in each pixel but I'm not sure how to do this in ModelBuilder. Do I have to somehow split each raster first, then add each one?

How would I iterate through each band of every raster?

Edit:Basically this is what I have so far. Right now it will take 3 bands of 1 raster and add the bands then go on to create a polygon of a certain reclassed area. I want to do this on a large number of rasters but am unsure of how to make it work.


texts: Mask image based on no-data pixels in another image with different projection - GDAL or Python

texts:
I've got a land cover image which covers the whole globe, and a satellite image which covers a very small part of the globe. These two images have different projections, co-ordinate systems and pixel sizes. Also, the satellite image is not fully 'square' within the image, it is on a slant, and there are significant no-data areas around the edge of the actual image.

I want to extract the pixels in the land cover image which intersect the satellite image.

All of this needs to be completely automatic - either via a Python script or a GDAL command-line. In the past when I've had a polygon outline of the image I've used a gdalwarp command like

gdalwarp -overwrite -of GTiff -cutline polygon.shp -crop_to_cutline landcover.tif output.tif

But I can't find a good automatic way to produce a polygon of the footprint of my satellite image, given than I don't want the footprint of the whole image, but just the part of the image that isn't No Data. So it doesn't look like that method will work. Theoretically I should be able to do a simple mask, but I can't work out how to do this when my images are of different resolutions.

Does anyone have any ideas how I should go about doing this?


texts: Trouble with Difference/Symmetrical Difference Tool - Empty output

texts:
I'm having issues trying to eliminate overlapping parts of two polygons.
The aim is to delete the parts of the pink layer that overlap with the green one. I tried using difference tool and symmetrical tool, however, both are producing an empty result (EPSG is the same)

Any ideas why this would be happening?
Also, if there are any alternative tools, can someone suggest them?

EDIT: I managed to fix it. There was a small intersection error in the layer that I checked with the Geometry Checker Plugin. After deleting that small error, the difference output was correct.


texts: Filling in holes in heightmap with lower resolution heightmaps?

texts:
I have this elevation map of the Verdal area in Norway, from the Norwegian mapping authority. The first one has a 1m resolution, but has some holes in the data. I'd like to fill in those holes with the corresponding area of the second height map, which has a 10m resolution, and merge them into one layer.

Data is from https://hoydedata.no/LaserInnsyn/ (Center coordinates ca 63.79142°N 11.43853°E, Nasjonal Høydemodell/National Elevation Model)

Does anyone know how to do this?

I am new to GIS software.


texts:
I have a vector layer containing the potential drainage lines that have been computed from a DEM. I want to be able to click on any point and identify all the lines that are connected to my point of interest. In effect, this means identifying all the streams that contribute water to the point.

I am sure this is not a novel problem. It would be similar to identifying all electric lines/houses that get power from a particular sub-station and other real life problems. I just can't seem to phrase my problem just right to get the correct response from Google.

I understand that a simple identification of all connected lines would not take care of upstream and downstream drainages. I guess that will have to figured out later.


texts: Creating Replica using ArcPY with check_out type

texts:
I am running the following command in python and it works fine no error but it does not generate a mdb in output path folder and give me the input path result as output child source

arcpy.CreateReplica_management("AHSAN.Point_1120","CHECK_OUT",
r"D:\Ahsan\January-16\Zaeem         \mdb_generation","CheckOut_Output_3","FULL","CHILD_DATA_SENDER","USE_DEFAULTS","DO_NOT_REUSE","DO_NOT_GET_RELATED","AHSAN.Point_1120","ARCHIVING")

<Result 'D:\\Ahsan\\January-16\\Zaeem\\mdb_generation'>

I need to create mdb in output path folder. I have tried the same thing with model builder but no difference.


texts:
I have a svg file:

https://drive.google.com/open?id=0B-u7dqmskzd5MEpuTmtZS3pDeHc

This is a single line that I want to convert into polygons.
If I understand correctly, I must convert it from svg to dxf first. Which I have done:

https://drive.google.com/open?id=0B-u7dqmskzd5eHZURVVsU2FwQU0

What next? I tried many solutions and nothing works.


texts: Excluding nodata value in band calculation with Rasterio?

texts:
How do I exclude no data values when doing a raster calculation with rasterio. I realise that the files are read into an array, and that "0" is a nodata value. However, when I do a simple sum like "image + 1" it applies the math to the nodata value as well. So what is the best way to 'ignore' nodata values in the calculation?

FYI: I am calibrating the Landsat bands to reflectance and this is what my current snippet of the band math looks like:

    if "B1.TIF" in band:
        print band
        image = "%s" % (band)
        with rasterio.open(image) as src:
            image_read = src.read(1)
        #np.seterr(divide='ignore', invalid='ignore')
        constant = 0.01745329251994444444444444444444 #Constant is calculated (3.14/180) which is converting the sun-angle to sun_radians which was suggested by WOlfgang
        set_mb = "{:.16f}".format(float(mult_band1))
        set_ab = "{:.16f}".format(float(add_band1))
        toa_1 = (float(set_mb) * image_read.astype(float)) + (float(set_ab))
        solar_z = np.cos((90-float(sun_elevation))*float(constant))
        toa_2 = (toa_1.astype(float) / solar_z) * 100000
        kwargs = src.meta
        with rasterio.open(outfile, 'w', **kwargs) as dst:
            dst.write_band(1, toa_2.astype(rasterio.uint16))
        #print image_read.mask, toa_2

Is there a simple way to ignore nodata values in the calculation?


texts: Libraries to render vector-based maps in the browser

texts:
Is there a open source JavaScript library like Leaflet or OpenLayers which can be used in a web or mobile application and loads vector-based maps instead of images? I only need the road information and some POIs.

The only library I found was polymaps - but still I would need a server or service which serves the vector SVG data (preferable from OpenStreetMap). Also there is the nice Android project called Mapsforge but I would need it for a web application or other mobile devices.


texts:
I have a map where each feature of a polygon shapefile will be filled with a different image. I have those images stored as PNG files. In my shapefile, each feature has a name which is the same as its corresponding png file.

There are over 180 features, and while I can select each feature separately and fill them with the raster image, I saw this function in the Advanced button and thought maybe this would be an easier way to do it. It requires an XML file.

But I have no idea what the structure of this xml file would be.

What are the tags?

How do you associate a file?

I cannot find any material on this anywhere online.


texts: Avoiding General function failure and deleting table field in ArcGIS Desktop?

texts:
I am trying to delete a table field: From ArcMap table view, I right click on the field heading and select Delete field. Unfortunately, I obtain the following message:

Failed to delete the selected fields.
General function failure [TABLE_NAME]

and the field remains.

How should I proceed to delete this field?


texts:
I have sets of lat/lon points for vortices in the ocean, and what I want to do is draw curves around them (they are roughly elliptical), and determine their eccentricity and orientation as they move through the ocean in time. I'm using python, and I'm on a mac.

I am having a really dreadful time navigating the sea (pun intended) of special libraries that are intended for this purpose.

Can anyone provide me with a simple list of instructions to do this?

I found some tools that point to Fiona and Shapely as easy-to-use tools, but I cannot for the life of me figure out how to get this array of lat-lon points into a .shp file which Fiona can open.

1) How do I create a .shp file out of these data points? If there is another filetype I should use, please let me know how. This is as far as I've gotten:

import Shapefile
w = shapefile.Writer(shapeType = 8) #8 is multipoint

From here I'm totally stuck, and can't find documentation that helps me past this. I don't need any additional features, literally just the points.

I have tried the following as well:
    w = shapefile.Writer()

for i in range(vortex.shape[0]):
    w.point(vortex[i, 0], vortex[i, 1])

for i in range(0, 327):
    print w.shapes()[i].points

w.save('sample_vortex.shp')

Vortex is a numpy array with dimensions (327, 2), with each row containing a lat and a lon value (duh).

This prints out a list of the correct points, so I'm on the right track. I found this snippet of code, which I think should plot my points:

import fiona
import shapely.geometry as geometry
input_shapefile = 'sample_vortex.shp'
shapefile = fiona.open(input_shapefile)
points = [geometry.shape(point['geometry'])
      for point in shapefile]

import pylab as pl
x = [p.coords.xy[0] for p in points]
y = [p.coords.xy[1] for p in points]
pl.figure(figsize=(10,10))
pl.plot(x,y,'o', color='#f16824')

But it actually does not, and yields a blank plot, which tell me, I think, that my .shp file isn't correct.

Any ideas?

2) Do you know how to draw these ellipses and get the eccentricity (and how much the major axis deviates from horizontal, while we're at it)?


texts: Segmenting rasters into sub-regions and sample using R without using extract()

texts:
Is there a more efficient way to sample from a sub-region (defined by a shapefile read into a spatial object) of multiple rasters without having to use extract() for each region on every raster?

Data: I have a large study area represented by approximately 50 raster layers, one per variable (e.g. elevation, slope, vegetation, etc...) my study area is divided into 10 sub-region

Methods: I need to draw random samples (using sampleRandom() or sample() depending on object) from each region of each raster; therefore approximately 500 samples.

Issue: I can do this relatively easy by using the extract() function for small datasets, but it will take forever with these data sets.  I am hoping there is a way to sample from within the area of the sub-regions without having to use extract()

Below is code (some thankfully borrowed from J. Evans) that shows the basics of what I am doing.
This works fine for small raster data sets, but trying to use extract() on large raster data sets fails.

require(raster)

r1 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r1[] <- runif(ncell(r1),0,100)
r2 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r2[] <- runif(ncell(r1),100,200)
r3 <- raster(ncol=10, nrow=10, xmn=0, ymn =0, xmx=100, ymx=100)
r3[] <- runif(ncell(r1),200,300)

region1 <- rbind(c(0,0), c(50,0), c(50,50), c(0,50), c(0,0))
region2 <- rbind(c(50,0), c(100,0), c(100,50), c(50,50), c(50,0))
region3 <- rbind(c(0,50), c(50,50), c(50,100), c(0,100), c(0,50))
region4 <- rbind(c(50,50), c(100,50), c(100,100), c(50,100), c(50,50))

polys <- SpatialPolygons(list(Polygons(list(Polygon(region1)), "region1"),Polygons(list(Polygon(region2)), "region2"), Polygons(list(Polygon(region3)), "region3"), Polygons(list(Polygon(region4)), "region4")))

plot(r1)
plot(polys, add=TRUE)

r_stack <- stack(r1, r2, r3)
names(r_stack) <- c("r1", "r2", "r3")

for(i in 1:length(polys)){
    region <- polys[i]
    for(k in 1:length(names(r_stack))){
        raster_region <- extract(r_stack[[k]], region)
        sample_values <- sample(unlist(raster_region), 10)
        print(sample_values)
        ### add samples to some matrix...
    }
}
### Do analysis on values...


texts: Create button to cycle through different markers in a GeoJSON file on click

texts:
I am trying to create an interactive Leaflet map of a hiking trail. Along this trail, I have a number of 'waypoints', represented as markers in a GeoJSON file:

var waypoints = [
   { "Name": "Cabane de Susanfe",
     "Altitude": 2095,
     "Text": "This is the first waypoint",
     "coordinates": [ 6.91246, 46.14939 ]
   },
   { "Name": "Rhône Valley",
     "Altitude": 457,
     "Text": "This is the second waypoint",
     "coordinates": [ 7.0458, 46.13417 ]
   },
   { "Name": "Cabane de Sorniot",
     "Altitude": 2064,
     "Text": "This is the third waypoint",
     "coordinates": [ 7.09167, 46.16176 ]
   },
   { "Name": "Cabane du Demècre",
     "Altitude": 2361,
     "Text": "This is the fourth waypoint",
     "coordinates": [ 7.07314, 46.18744 ]
   },
   ]

I would like to create a button in a separate div which would cycle through these waypoints and display information about each in a div on the side, like a map slideshow. So far I have managed to create a button that cycles through the elements in the 'waypoints' array, but how do I get the map to pan to the new marker and display information? Here is the jquery code I have so far to iterate through the array elements:

$(document).ready(function(){
var i = 0;

$("#button-next").on('click', function(){
     i = (i + 1) % waypoints.length;

     console.log(waypoints[i]);
});
});

The effect I am looking for is similar to CartoDB's odyssey.js slideshow or the jQuery slidemapper plugin, but I figure that there is probably an easier solution than completely rebuilding the map using one of these libraries.

I am a beginner with leaflet and javascript in general.


texts:
Using QGIS 1.8, I've got a shapefile with an integer field called "AREA" and I'd like to sum the values which are in this shapefile.
Does anybody know how to?