texts: Adding ArcGIS JSAPI code to Existing Angular app?

texts:
I am new to web dev. coding.  I managed to scrape together some code for a simple web GIS app. using ArcGIS Javascript API version 3.xx. The ArcGIS JSAPI app has some simple functionality like a click event that return some value back from the layers.  There is an existing Angular(version 7) app running in our environment that I inherited from my other colleagues, with different modules inside of the src-->app folder already.

Can I just simply create another module inside the src-->app folder of the Angular app, maybe called it "gis-map", and copy and paste the html/js code from the ArcGIS JSAPI app into the gis-map.component.html, " ".ts, " ".css?

Will this approach work?

And could I expand that later on, to add a Angular click event on the map layers and get back the result?


texts: Reading data from infowindow in Cartoons?

texts:
Is there a way to read data from a CartoDB infowindow and write it, e.g. into a form?

For example: In a city plan with parcels, the user clicks on a parcel and the infowindow opens with information on this parcel. Now I think about a “Fill in” button that allows writing the attributes of the parcel selected into a form.

Any hints?


texts: How to select features from different WFS layers?

texts:
I'm trying to select features from different WFS layers but unfortunatelly it doesn't works in my code.

This is my code:

var select = new OpenLayers.Layer.Vector("Selection", {styleMap: new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])});
var hover = new OpenLayers.Layer.Vector("Hover");

map.addLayers([hover,select]);

controles = new OpenLayers.Control.GetFeature({
    protocol: new OpenLayers.Protocol.WFS({
        url: "http://demo.opengeo.org/geoserver/wfs",
        featureType: ["states","world"]}),
    box: true,
    hover: true,
    multipleKey: "shiftKey",
    toggleKey: "ctrlKey",
});

controles.events.register("featureselected", this, function(e) {
      select.addFeatures([e.feature]);
});

controles.events.register("featureunselected", this, function(e) {
      select.removeFeatures([e.feature]);
});

controles.events.register("hoverfeature", this, function(e) {
      hover.addFeatures([e.feature]);
});

controles.events.register("outfeature", this, function(e) {
      hover.removeFeatures([e.feature]);
});

map.addControl(controles);
controles.activate();

As you can see, I've put in the featureType two layers but only the first one can be selected. What am I doing wrong?


texts: How do I import multiple tif images and save it as one tiff file with skimage or gdal in python?

texts:
I am new to Python and I am trying to import multiple tif images so that I can do image processing with skimage. I want to do image processing on a Landsat image to define linear features for a project that has to be in before the end of the year.

I am using a Landsat 8 image that consists of 11 bands (individual tif files), a bqa file and a metadata.txt file.

I can import a band one by one, but I don't know how to stack the bands into one tiff file that also contains the metadata.txt file.

This is what I have done so far:

import os filename = os.path.join('path', 'band1.TIF')
import skimage from skimage
import io bandone = io.imread(filename)

So my image has been successfully imported as converted to a nparray.
I have to repeat this to import all the bands.

What is the best way to save these bands together?

skimage.external.tifffile.imsave('file, data, **kwargs) ?
skimage.io.imsave(fname, arr, plugin=None, **plugin_args) ?
TiffWriter ?

Or should I do something else completely?

Please elaborate on the details (such as how to define **kwargs for instance).


texts: FitBounds after L.GeoJSON.AJAX has loaded

texts:
I am making a Leaflet 1.3.4 map that will show members based on category. I have a php file "cs.php" that queries a MySQL and returns an option menu of distinct categories. When a category is selected it makes the selected value available to the "showUser" function. That function should clearLayers() then then bring in the appropriate GeoJSON data displaying markers with popups. All of that works fine, my problem is trying to get "fitBounds" to to zoom to fit the markers on the map. I have tried multiple scripts. Whenever fitBounds is part of the showUser function the map does not show and not errors are in the console. I don't understand what I have done wrong or forgotten or something.

UPDATE: I updated the code with suggestion to use "geojsonLayer" by @Marc_Pfister suggested in the .on data:loaded and that fixed it. I also changed to data:progress per @Stefan_Fairphone.

Here is the working example:

Here is an example of the GeoJSON that "GeoJSON.php" is outputting:

{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"Point","coordinates":[-94.905793,39.416633]},"properties":{"Member_Name":"Angels Rock Bar","Contact":"Pedro Lopez","Address":"1323 Walnut","City":"Kansas City","Phone":"(816) 896-3943","Email":"info@angelsrockbarkc.com","Website":"www.angelsrockbarkc.com"}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-94.541461,39.051003]},"properties":{"Member_Name":"Bar K","Contact":"Leib Dodell","Address":"501 Berkley Parkway","City":"Kansas City","Phone":"(816) 225-6961","Email":"wag@barkdogpbar.com","Website":"www.barkdogbar.com"}},{"type":"Feature","geometry":{"type":"Point","coordinates":[-94.583506,39.098954]},"properties":{"Member_Name":"The Indie on Main","Contact":"Jamie Pitts","Address":"1228 Main St.","City":"Kansas City","Phone":"(816) 283-9900","Email":"","Website":""}}


texts: how to set default style from particular workspace using gsconfig?

texts:
I'm able to set default style for any layer using below code

cat = Catalog("http://localhost:8080/geoserver/rest", "admin", "geoserver")
layer = cat.get_layer("layer_name")
layer._set_default_style("style_name")
cat.save(layer)

Now I want to set default style from particula workspace (say states),How I can able to do?
I tried as

cat = Catalog("http://localhost:8080/geoserver/rest", "admin", "geoserver")
layer = cat.get_layer("layer_name")
layer._set_default_style("states","style_name")
cat.save(layer)

The error is: TypeError:_set_default_style() takes exactly 2 arguments (3 given).


texts: How to select more than one feature using values?

texts:
What should I do if I want to select more than one feature using values ?

For selecting a single feature, I enter an expression like this which works:

"LINK" = 28

Assume LINK as a field in the above example. Now what if I want to select more fields like that. I expect something like this:

"LINK" = [28, 56]

The above expression doesn't work (obviously). Is there a way to achieve this in Quantum GIS ?


texts: Automate a process in QGIS using graphical modeler

texts:
I have tasks that I need to automate:

Dissolve line features by attribute - (Dissolve tool)
Split the dissolved layer into single parts - (Multipart to singleparts tool)
Extract different layers based on an attribute(values are integers e.g. from 1 to 6) - (Split vector layer tool) - here the result is many layers depending on unique attributes. The results are saved in a directory(folder).
Intersect one extracted layer with another layer - (Split line with lines). Here only two layers can be intersected. there is a pattern to intersect for example where attribute = 1, then intersect with layer whose attribute =2. Therefore only an increment of 1
Merge the resulting layers.

As input I have a stream network that is quite segmented. This process ensures that I have the correct number of segments based on my requirements.

I am able to set up the tools until number 3. I have set up the input parameters as the layer and the field to dissolve. The result of number 3 is many different layers that are stored in an output directory. I have no clue on how to access the output directory and set up the Intersection (number 4).


texts:
I am working on a script that will populate 3 letters into a field. I have to use letters due to field limitations and to make sure that I have enough combinations. I'm not really sure where I am going wrong but I don't get the repeated values filtered out consistently. It seems to happen most of the time, just not all. Any one see anything wrong in my code?

import random, string, arcgisscripting, time

gp = arcgisscripting.create(9.3)
gp.OverwriteOutput = 1

# Lets set up some variables
layer = "Layer Goes Here"
field = "Code"
length = 3 # Declares the number of character to be created.

uniquelist = [] # Empty list to store already created characters.

##########Define the random character generator##############

def ranchar(length):
    myrg=random.SystemRandom

    characters=["a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w","x","y","z"]
    pw=string.join(myrg(random).sample(characters,length))
    randomval = (pw.replace(' ', '')) # Removes some spaces between the letters.

    return randomval

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

cur = gp.UpdateCursor(layer)
row = cur.Next()

while row:
    randomval = ranchar(length) # Grabs a random character.

    if randomval in uniquelist: # check to see if it already exists.
        print str(randomval) + " already exists, generating new value"
    else:
        uniquelist.append(randomval) # If the value isn't in the list, adds it.
        print uniquelist[-1] # Just for display.
        row.SetValue(field,uniquelist[-1]) # Grabs the last value from the list and writes it.
    cur.UpdateRow(row)
    row = cur.Next()

# Delete cursor and row objects to remove locks on the data
#
del cur, row

Thanks,

-Mike


texts: Cleaning up GPS coordinates with Google imagery?

texts:
I have a set of noisy GPS coordinates that I would like to clean up using Google Earth/Street View. Examples are shown below. Basically, I would like to be able to click and drag my existing coordinates (with all of their attributes) to what I think is the correct location based on what I can see in the imagery.
Is there a standard/easy workflow for doing this? For example via QGIS or Google Earth engine? Currently I am using Google maps, clicking on the desired location, and copying the new GPS coordinates manually into a spreadsheet. But this seems very inefficient, I would prefer to just click the coordinates and drag.
Thanks!


texts: In PostGIS, why does ST_Covers return f on this query?

texts:
In general this query functions fine. In my PHP application I tested it with a large number of points, and they all seemed to be true or false as appropriate (when I pasted data into two online WKT converters). But there's a significant number of false negatives. Here's the query using a point PostGIS thinks doesn't go in the rectangle even though I'm quite sure it should:

SELECT
    ST_Covers(geog_poly, geog_pt)
        As poly_covers_pt
            FROM (
                SELECT
                    ST_GeogFromText('POLYGON((-53 83,-53 42,-141 42,-141 83,-53 83))')
                        As geog_poly,
                    ST_GeogFromText('POINT (-100 48)')
                        As geog_pt
                )
                    AS result;

That returns poly_covers_pt as f, even though 42 < 48 < 83 and -141 < -100 < -53. Also, mapping the point and the region with online tools shows the region to contain all of North Dakota and the point to be in North Dakota.

By contrast, I can change the 48 to a 60 and this returns a t, as it should:

SELECT
    ST_Covers(geog_poly, geog_pt)
        As poly_covers_pt
    FROM (
        SELECT
            ST_GeogFromText('POLYGON((-53 83,-53 42,-141 42,-141 83,-53 83))')
                As geog_poly,
            ST_GeogFromText('POINT (-100 60)')
                As geog_pt
        )
            AS result;

What's going on?


texts: Shapefile location change repair. HResult E_FAIL

texts:
I'm working on a ArcMap plugin using ArcObjects that will change the location of a shapefile that is broken. However, I'm getting HResult Error E_Fail whenever I try opening the feataure class. This is the specific HResult value: -2147467259

Here's the setup. I have a MXD file with a shapefile that I created. This shapefile's location has been moved from the MXD file expects it to be. It expects it to be in C:\usr\Desktop, but I have moved it to C:\usr\Desktop\tmp. The shapefile is called "New_Shapefile.shp"

The goal is to repair the location of the shapefile. Here is the current code I have for shapefile replace. The following fails whenever I try opening the feature dataset.

IFeatureLayer featureLayer = (IFeatureLayer)pLayer //pLayer is the current ILayer;
IFeatureClass oldFeatureClass = featureLayer.FeatureClass;

IPropertySet propertySet = new PropertySetClass();
propertySet.SetProperty("DATABASE", "C:\usr\Desktop\tmp");
IWorkspaceFactory workspaceFactory = new ShapefileWorkspaceFactoryClass();

try
{
    IWorkspace workspace = workspaceFactory.Open(propertySet, 0);
    //I've also tried the following for opening up a workspace:
    //IWorkspace workspace = workspaceFactory.OpenFromFile("C:\usr\Desktop\tmp", 0);
    IFeatureWorkspace fWork = (IFeatureWorkspace)workspace;

    //program hits exception on line below
    IFeatureDataset newFeatureClass = fWork.OpenFeatureDataset("New_Shapefile");

    IMapAdmin2 mapAdmin2 = (IMapAdmin2)map;
    featureLayer.FeatureClass = newFeatureClass;
    mapAdmin2.FireChangeFeatureClass(oldFeatureClass, newFeatureClass);

    ESRI.ArcGIS.ArcMapUI.IContentsView contentsView = ArcMap.Document.CurrentContentsView;
    contentsView.Refresh(null);
}
catch (Exception e)
{
    ErrorMessage("uh oh something went wrong :(", e);
}

This code was based off of: http://edndoc.esri.com/arcobjects/9.2/ComponentHelp/esriGeodatabase/IFeatureWorkspace_OpenFeatureClass_Example.htm

Edit: This is for ArcMap 10.2

Edit: This code block works, I just had to rebuild the solution and it's working now. I'm not sure why it required a rebuild though. I've tested this for geodatabase files as well, just change the ShapeFileWorkspaceFactoryClass() to FileGDBWorkspaceFactoryClass()


texts: Overpass query: 2 almost identical search queries - one works, one not. Why?

texts:
I'm using a tool called "Wikivoyage Districitifier" to visualize OSM boundary relations. The query below is supposed to get all relations with certain admin_level values within a defined bounding box from OSM.

For Sydney it works perfectly okay (Overpass turbo):

[out:json][timeout:60];(relation[admin_level](if:(t["admin_level"] >= 5 &&t["admin_level"] <= 11))["wikidata"](-33.86400143885985,151.2000692923067,-33.85388041541516,151.21772895646075););out body;>;out skel qt;

For Auckland, the following query returns an empty result, which Wikivoyage Districitifier obviously is unable to process, as it only reports an error message: "Bounds are not valid". Overpass turbo):

[out:json][timeout:60];(relation[admin_level](if:(t["admin_level"] >= 5 &&t["admin_level"] <= 11))["wikidata"](-37.025032151631635,174.4815813201933,-36.71301768775455,175.04669057312287););out body;>;out skel qt;

Does anybody know what I do wrong?


texts: Problem with grid labels in grid composer

texts:
I am having some trouble with the print composer. First I set a fixed scale the proper scale for my map. Second, I set the width and height. Third I set the fixed X and Y minimum for my map.

For instance, I have a map that I want at a 1: 5,000 scale. Width is 240mm which would display 1,200m. Height is 160mm which would display 800m. X minimum is 481,000 and Y minimum is 2,136,000.

When I add a grid, I expected to see my minimum X and Y in the lower right corner and on the Top right corner. Nevertheless, the labels do not appear. I can only see labels for 482,200 in the Lower and Top corners, no label for my Y minimum in the lower left corner but the 2,136,800 at the top right corner.
I wonder if I am missing a step or some option.


texts: Dotspatial.Projections with custom Proj.4 string

texts:
I have lat,long values that I need to convert to x,y in a custom Lambert Conformal Conic projection defined as +units=m +a=6371229 +b=6371229 +x_0=0 +y_0=0 +proj=lcc +lon_0=0 +lat_2=25.0 +lat_1=25.0 +lat_0=25.0. I'm trying to use Dotspatial.Projections library in C# to basically replicate a process I have working in Python using Pyproj. Unfortunately, my C# function produces 0,0 which is not correct. Please find below the code:

string pj_grid_str = @"+units=m +a=6371229 +b=6371229 +proj=lcc +x_0=0 +y_0=0 +lon_0=0 +lat_2=25.0 +lat_1=25.0 +lat_0=25.0";
string pj_latlon_str = @"+proj=latlong +a=6371229 +b=6371229";
var pj_grid = DotSpatial.Projections.ProjectionInfo.FromProj4String(pj_grid_str);
var pj_latlon = DotSpatial.Projections.ProjectionInfo.FromProj4String(pj_latlon_str);

double[] ll = { -114.0, 51.0 };
double[] z = { 0.0 };

DotSpatial.Projections.Reproject.ReprojectPoints(ll, z, pj_latlon, pj_grid, 0, 1);

 System.Console.WriteLine(ll[0].ToString() + "," + ll[1].ToString());   //produces "0,0" no matter what input I provide
// should produce "-1489098.19006,3107428.30237"

Again, using a similar procedure in Pyproj with the same proj.4 strings works just fine and produces correct results.

P.S. I seem to have found a way to do this using Proj.NET library, but am still puzzled as to why Dotspatial can't do it.


texts: refresh programmatically created composer

texts:
I am trying to create a composer using PyQGIS. Everything works as expected. The only difference compared to the manual creation of a composer is, that I am not able to update the layers (by pressing "Refresh View") in the Composer GUI after creating the composer by script which works fine when I create the composer the GUI-way.

My code so far:

from PyQt4.QtXml import QDomDocument

template_filePath = "PATH_TO_COMPOSER_TEMPLATE"
template_file = file(template_filePath)
template_content = template_file.read()
template_file.close()
document = QDomDocument()
document.setContent(template_content)
composition.loadFromTemplate(document)


texts: QGIS3 saved Raster layer is "skewed" or aspect ratio is wrong

texts:
I'm trying to export two layers based on the current view or "extent":

A black and white "raster" derived from some vector data
A map layer (Bing, Google, osm, etc)

The issue I'm having is that the raster layer (#1 above) has an aspect ratio that seems to be saved wrong, causing my data to not overlap correctly, and an obvious stretching (images below)

Warped Image

Expected Image

Save Dialog settings

I'm exporting the extent by right clicking on the layer, selecting Export > Save As, then setting the extent to the map view and saving as GeoTIFF. I'm saving it as the same format it came in as (WGS 84 EPSG:4326), and everything is lined up perfectly in the map view.

When using: Project > Import/Export > Export Map to image, it works perfectly :shrug:


texts: Error 000628: Cannot set input parameter maximum_pyramid_levels

texts:
I'm trying to execute a script that will create a mosaic dataset that includes Landsat imagery and I keep running into this error.  I'm fairly new at writing Python scripts, so I'm sure I'm just missing something small, but I honestly cannot figure this one out after searching around a lot for a fix.  My code is below:

# Imports
import arcpy
import os
from arcpy import env
env.overwriteOutput = True

# Define Variables
image_dir = r"E:\GIS 4090\Week 11\In Class Exercises\Data\Landsat_p114r75"
output_gdb = r"E:\GIS 4090\Week 11\In Class Exercises\Results\Challenge.gdb"
mosaic_name = "Multispectral"

arcpy.env.workspace = image_dir
image_list = arcpy.ListRasters()
print(image_list)

# Create File Geodatabase
# Use os to get the output_gdb+dir ad output_gdb_name from output_gdb.
output_gdb_dir = os.path.dirname(output_gdb)
output_gdb_name = os.path.basename(output_gdb)
if arcpy.Exists(output_gdb) != True:
    arcpy.CreateFileGDB_management(output_gdb_dir, output_gdb_name)
else:
    print("File Geodatabase exists!")

# Create Mosaic Dataset
# Tool help at http://pro.arcgis.com/en/pro-app/tool-reference/data-management/create-mosaic-dataset.htm
projection = arcpy.SpatialReference(32750)
arcpy.CreateMosaicDataset_management(output_gdb, mosaic_name, projection)
print("Mosaic Dataset created.")

# Add rasters to mosaic dataset
# Use os.path.join to create the "path_to_mosaic" variable
# http://pro.arcgis.com/en/pro-app/tool-reference/data-management/add-rasters-to-mosaic-dataset.htm
rasttype = "Landsat 7 ETM+"
updatecs = "UPDATE_CELL_SIZES"
updatebnd = "UPDATE_BOUNDARY"
updateovr = "UPDATE_OVERVIEWS"
buildpyr = "BUILD_PYRAMIDS"
calcstats = "CALCULATE_STATISTICS"
buildthumbs = "BUILD_THUMBNAILS"
opdesc = "6 band Landsat imagery for analysis"
path_to_mosaic = os.path.join(output_gdb, mosaic_name)
arcpy.AddRastersToMosaicDataset_management(path_to_mosaic, rasttype, image_list, \
updatecs, updatebnd, updateovr, buildpyr, calcstats, buildthumbs, opdesc)
print("Rasters added to mosaic.")

# Not sure if the following code section is needed:

#Calculate Statistics on Mosaic Dataset
skip_factor_x = "5"
skip_factor_y = "5"
arcpy.CalculateStatistics_management(os.path.join(output_gdb, mosaic_name), \
skip_factor_x, skip_factor_y)
print("Statistics Calculated.")

print("Done.")


texts: How do you import a .SHP into Sketch-Up?

texts:
Does anyone know how to use a .SHP in Trimble Sketch-Up?

I have tried exporting to KML and CAD from QGIS, but the KML file doesn't read in Sketch-Up.


texts: Noding a network with grade separation - open source solution?

texts:
Various network data standards (such as those used by Open Street Map, pgrouting) employ "nodes for junctions, intersections for bridges/tunnels" rule.  That is to say, links are considered disconnected if they merely intersect but joined if they have coincident endpoints.

Other standards use a rule of "nodes wherever links cross", even if there is no node in practice (e.g. bridge, tunnel).  The ambiguity thus caused is resolved by grade separation data.

Is there an open source tool for converting the former into the latter?


texts: How do I reformat numbers to only show a certain number of significant digits in CartoDB?

texts:
I need to round numbers that are in my table in CartoDB, so that in the Map visualization infowindows I only show results to two decimal digits rather than the long decimals that are currently showing.

I tried using SELECT round(column,10,2) FROM table as suggested on this webpage How can I show numbers as percentages in CartoDB?

but the SQL query results in no records, and I get the following error:

function round(double precision, integer) does not exist

Is there another way to reformat the numbers so they display fewer decimals in the CartoDB infowindows without editing the source data?


texts: How to use Zoom Slider like the older version in Arcgis Javacript API 3.3?

texts:
I'm using Arcgis Javacript API 3.3. In this version zoom slider look like this.

And here, this is the older version.

Is there any way to make it like the older version?


texts: What criteria should be used to select a project CRS in QGIS?

texts:
I've been using QGIS for a little while now, but I am still a bit confused about coordinate systems. There is plenty of guidance about how to change project and layer CRS, on-the-fly transformation etc, but I can't find anything which explains fairly simply why you should use a particular CRS in preference to another. Does it depend on the data sources (if any) which I intend to use? Or the size of the area to be covered (eg a few square kilometres vs a whole country)? Or the locale?

I've tried searching for tutorials etc which can shed a bit more light, but haven't come up with anything yet.


texts: How to keep an OSM routing topology updated?

texts: I have a South America import and an osm2po topology table. My database is updated daily using osmosis but how to keep the topology table updated without make all import over again?


texts: Importing shapefile into GRASS, and multiple new rows appear in attribute table

texts:
I am importing a GADM Global Administrative Areas shapefile into GRASS (using QGIS 2.8.2 interface).

When I do this, multiple additional rows are created in the attribute table.

It goes from this:

To this:

This makes it difficult to work with. I am using v.in.ogr with and without the "-c" flag, but this doesn't make a difference.

Is there anything I can do to keep the attribute table in the original form?

There is a discussion on importing GADM into GRASS stating that importing causes it to crash, but it doesn't help here. Further I am importing a very small subset of GADM (~20 polygons, some with edges touching). There is documentation on simplifying GADM, but this is beyond my capabilities.


texts: How to make two feature layers congruent using ArcGIS for Desktop?

texts:
I have two feature layers (polyline/polygon) which should be congruent but there's an offset of about 220 meters.

The Georeferencing tool only works with raster layers.

How do I do the equivalent with a feature layer?


texts: Computation Timed out error when creating a cloud-free composite of Sentinel-2

texts:
I am trying to create a cloud free-composite that covers my ENTIRE study area using Sentinel-2 data. But, everytime I run this code, I get a Computation Timed out error. I am simply loading an image which is a composite of 10 images or lesser.

Link to code:(features can be accessed easily) https://code.earthengine.google.com/0d9ffb91268960807edc86bec867ee65

And code below

    //Set dates
var year = 2015;
var start_date = '-07-01';
var end_date = '-08-30';
var end_year = year+1;
print (year.toString() + start_date,end_year.toString() + end_date);

// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = ee.Number(2).pow(10).int();
print (cloudBitMask); //1024
var cirrusBitMask = ee.Number(2).pow(11).int();
print (cirrusBitMask); //2048

var cloudmask = function (image) {
  var qa = image.select('QA60');
 // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask);
};

// Create clipping function
var clipper = function(image){
  return image.clip(WG);
};

//// Application
var clipped = sentinel2.map(clipper);
var filtered = clipped.filterDate('2015-07-01','2018-08-30')
                      .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE',5))
                      .map(cloudmask)
                      .sort('system:time_start',false)
                      .limit(10);

// Create RGB composite. I use the median of the 10 most recent images
// to enhance the changes that really no clouds are present and the
// complete ROI is covered.
var median = filtered.reduce(ee.Reducer.median())
                   .select(['B4_median','B3_median','B2_median']);
Map.addLayer(median);

Any suggestions / solutions?

Before you berate my post, YES, I did take a look at the solution here and it did not help - How do I create a Sentinel-2 cloud free, cloud-shadow free composite or scene on Earth Engine (GEE)?


texts: Using "Difference" geoprocessing tool on Big Data taking a very long time

texts:
I am using QGIS Bonn 3.2.3.

I've recently been using the Difference function to calculate a list of points and it is taking a terribly long time. I was wondering if you could point out any obvious mistakes or better ways of solving my problem.

My end goal is to create a list of properties in Scotland that are likely to have a poor internet connection. I do it by considering two sets of points: Properties and Cabinets. There are about 3.5mln properties and 20k cabinets. The data sources are a .csv text file

Then I draw a 1km buffer around the cabinets and run there difference function to find the properties which are more than 1km away from cabinets.

Sounds pretty simple and everything works except that the calculation took 25 hours to finish and QGIS wasn't even using 10% of my computer's resources. And I am finding myself having to repeat a similar calculation for telephone exchanges and cabinets as well as for different distances in other UK countries. Current overall computation time estimate is over a week. And that's assuming I get it right the first time

My question: Is there a way of making QGIS use more of my computer's resources to finish the calculation faster or is there some other functionality I might use? Please share some other tips and tricks that could help me out.

At the moment, I feel like QGIS isn't the right tool for my task. I am considering writing a python script to compute this problem but (if possible) I would much rather have it run in QGIS.


texts: Cannot create postgis extension to postgres database even though postgis is installed

texts:
What I do:

Windows 7, Pgadmin, I select my database, go to extensions, select "New Extension", select "postgis" from List. Then, I get: "An error has appeared: Error: could not load library: C:/Program Files/PostgreSQL/9.4./lib/postgis-2.1.dll". However, when I look into this folder, the file definitely exists.

What can be the problem?


texts: How to replace raster values in one Raster with another

texts:
I have two raster files in ArcGIS 10.1, raster 1 and raster 2 all with values and similar spatial resolutions, and overlay in some areas.

I want to replace values in raster 1 with those in raster 2 where there is an overlay.


texts: Mapping coordinates from one to another location?

texts:
I have some MapInfo files in my workspace in FME.

The coordinate system of my data is "LL-GDA94", this coordinate system is related to Australia.

My problem is the coordinates of all data in TAB files show somewhere in Japan (but not in Australia).

Is there any solution that can map the coordinates in Japan to Australia?

For example in the following image you can see the coordinate system in _coordsys column and the coordinate of a point in the red circle :

33.922691 (X), 136.059935 (Y)

If I copy that coordinate to the google map, it's somewhere in Japan:

Where it should be -33.830366 (S), 151.22157 (E)


texts: How to convert a GRIB file to CSV/ASCII?

texts:
I have downloaded the ERA5-Land hourly-ECMWF climate reanalysis data from Copernicus. The data come in .grib format. It would be easy for me to manipulate the data in .csv or ASCII format. The idea is to have data presented in a "tabular" form with the columns: time/longitude/latitude/precipitation/temperature.
I have done some reading and I have found the following solutions that I cannot implement.

wgrib2 : I can't seem to find the place to download it (or link not working)
Pygrib : I cannot use this python library because I am using windows 10
google earth engine: I cannot find a script available online.

Using windows 10 as my OS:

What python library (with script if possible) can I use to convert grib to csv or ASCII?
Or where can I find a script to do the conversion on google earth engine?
Or any software that can allow me to do such convertion?


texts:
I want to create map overlays for a TomTom GPS. Relief shading, ortophoto etc as map background while navigating in Map View mode.

Instructions for creating custom overlay can be found here

So far I have download relief shading and ortophoto imagery from wms server for whole Estonia. (Relief shading and ortophoto would be separate overlays, I don't want to combine them.) That's over 12000 .jpg files for each. Files are 2000x2000 pixels (10x10km) in size (except tile files from last row and column, which are smaller).

What I need is to create one .sat file per each .jpg file.
.sat file is simple text file, similar to world file used by ESRI.
.sat file has same name as the .jpg it reffers to (except the extension)
.sat file contains 9 rows of information:

Line 1 - filename of the image file.
Line 2 - GPS coordinate of the top left corner of the image (longitude).
Line 3 - GPS coordinate of the top left corner of the image (latitude).
Line 4 - GPS coordinate of the bottom right corner of the image (longitude).
Line 5 - GPS coordinate of the bottom right corner of the image (latitude).
Line 6 - Minimum zoom level for image to be visible (min = 0).
Line 7 - Maximum zoom level for image to be visible (max = 65,535).
Line 8 - Width of image file in pixels.
Line 9 - Height of image file in pixels.

So how could I generate such files automatically? I'm using ArcGIS.


texts: Offset Overlapping Routes w/out Schematics

texts:
I'm currently working on showing route frequency and productivity for a transit system. The data is by route and I'd like to have the routes offset to better visualize the data.

I've researched this and have found that I can use a combinations of geometric networks and the ArcGIS Schematic extension. Is is possible to do this without the Schematic extension?

I know I could simply change the symbology of every route line, but with them routes and going in multiple directions, simply changing the offset in the line's symbology really wouldn't work.

In the past I've manually offset each route and I'd really like to not have to do that again.


texts: Calculating field AVG from B layer 3 nearest points to each A layer point

texts:
I wonder how to calculate statistics (in my case average) of one point layer field (let say B layer) and assign them to each point of A layer.
Statistics should be calculated from 3 nearest B layer points that are closest points to the point from A layer.
My software is QGIS 3.4.15.

For the above example:

Let say that I want to calculate avg of 'Price' attribute that contains B layer.
The AVG(B2.Price, B3.Price, B4.Price) should be assigned to point A1 of A layer.
The AVG(B5.Price, B6.Price, B7.Price) should be assigned to point A2 of A layer.
The AVG(B5.Price, B8.Price, B9.Price)  should be assigned to point A3 of A layer.
It would be great to assign also the distance to the 3rd point.


texts: Rule based display labels, show different value with first is null

texts:
I am attempting to create a mapping system for navigation at my small fire department using QGIS and a GPS.  I am attempting to show spaces in a trailer park or apartment building instead of the street address.  So instead of the label 1234, 1234, 1234, 1234, over and over, it would read (space: 1, 2, 3, 4, 5) Some of these are very close together so I was hoping that if value "space Number" is present then value "address" is absent.  Also, if there is no space number it would just read the address.

I attempted the "ADDRESS || '\n' || SPACE" and when I do so it only shows the labels that have both values.  I just want one or the other, how could I write that expression.


texts: raster2pgsql with the -R flag (remote dataset) results in invalid table

texts:
I'm trying to setup a PostGIS table with NED raster tiles stored remotely. My raster2pgsql command works fine and results in a valid table without the -R flag, but when I use the -R flag (for an 'out-of-db' raster), the resulting table does not have a correctly defined CRS and the postGIS table (raster layer) is not visible when loading it into QGIS. This issue with the -R flag happens to both local dataset, as well as a COG located on an S3 bucket.
Is it possible to create a valid PostGIS raster table with just references (remote raster datasets)? Or am I going about this all wrong? This was also an interesting project that seems related, but i'm not sure where to take it.
Step 1: I have two NED 1" tiles that I'm testing and have a .txt doc listing these file locations:
Local (local_ned_1.txt):
NED_1/float/n48w114/floatn48w114_1.flt
NED_1/float/n47w114/floatn47w114_1.flt

S3 (s3_ned_1.txt):
/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n47w114/USGS_1_n47w114.tif
/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n48w114/USGS_1_n48w114.tif

Step 2: build a .vrt from the .txt files. This resulting .vrt works fine in QGIS for both local as well as the S3 cogtiff source:
Local: gdalbuildvrt -input_file_list local_ned_1.txt local_ned_1.vrt
S3: gdalbuildvrt -input_file_list s3_ned_1.txt s3_ned_1.vrt
Step 3: Analyze .vrt contents:
Local (just a sample showing paths):
<ComplexSource>
      <SourceFilename relativeToVRT="0">/Volumes/RASTER/NED/NED_1/float/n48w114/floatn48w114_1.flt</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="3612" RasterYSize="3612" DataType="Float32" BlockXSize="3612" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="3612" ySize="3612" />
      <DstRect xOff="0" yOff="0" xSize="3612" ySize="3612" />
      <NODATA>-9999</NODATA>
    </ComplexSource>
    <ComplexSource>
      <SourceFilename relativeToVRT="0">/Volumes/RASTER/NED/NED_1/float/n47w114/floatn47w114_1.flt</SourceFilename>

S3 (just a sample showing paths):
<ComplexSource>
  <SourceFilename relativeToVRT="0">/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n47w114/USGS_1_n47w114.tif</SourceFilename>
  <SourceBand>1</SourceBand>
  <SourceProperties RasterXSize="3612" RasterYSize="3612" DataType="Float32" BlockXSize="256" BlockYSize="256" />
  <SrcRect xOff="0" yOff="0" xSize="3612" ySize="3612" />
  <DstRect xOff="0" yOff="3600" xSize="3612" ySize="3612" />
  <NODATA>-999999</NODATA>
</ComplexSource>
<ComplexSource>
  <SourceFilename relativeToVRT="0">/vsicurl/https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/n48w114/USGS_1_n48w114.tif</SourceFilename>

Step 4: Load .vrt into postGIS with raster2pgsql:
Use the -R flag, this results in a table that does not show the raster layer when the table is loaded into QGIS, and has an invalid CRS:
raster2pgsql -c -s 4269 -R -I -Y -C /Users/username/local_ned_1.vrt public.ned1_local | psql -h localhost -p 5432 -U postgres -d database_name

DON'T use the -R flag, this results in a table that loads the actual raster data into postGIS (not what i'm after) but loads fine into QGIS and has a valid CRS (EPSG:4269)
raster2pgsql -c -s 4269 -F -I -Y -C /Users/username/local_ned_1.vrt public.ned1_local | psql -h localhost -p 5432 -U postgres -d database_name


texts: Printing layer name using label element in Print Composer in QGIS 3.12?

texts:
I'm trying to print the layer name using a label in print composer in QGIS 3.12.

It's not responding to layer_property(@layer_id, 'name'), or [%layer_property( @layer_id, 'name')%] or [%@layer_name%]. Any suggestions?

I have many layers and just one print layout, it is updating data with current selected layer, but not the layer name, the layout is not printing it.


texts:
I'm using MOD09GQ product of MODIS to calculate vegetation indices in Erdas 2014, and I was able to extract NDVI using Erdas Model Builder.

What is the formula I should use in model builder for calculating Normalized Difference Water Index (NDWI) in Erdas? Do i need to download another products of MODIS for computing NDWI?


texts: check if the district belongs to the province

texts:
I'm building map viewer using d3js and data from WFS service.
In Poland we have 3 levels of administrative divisions:

province
district
community

Using WFS service I can get exact boundaries of each level.
For example for provinces I'm using http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getFeature&version=1.1.0&service=WFS&typename=gmgml:WOJEWODZTWA
For districts: http://sdi.geoportal.gov.pl/wfs_prg/wfservice.aspx?request=getFeature&version=1.1.0&service=WFS&typename=gmgml:POWIATY

Then using ogr2ogr I'm converting those GML files to TopoJSON and displaying using d3js. This part works fine.

Right now I need to know which district belongs to which province, but unfortunately this isn't available in file that comes from WFS.
I need this to be able to filter my map by provinces - display only districts for specific provinces.

I know that I can open SHP that file in QGIS and add attribute for every district, but I would have to add this to 2000+ items (before that checking manually if that district belongs to that province), and I would like to avoid that extra work.

Can this be done automatically? Using some kind of software?
I have only tiny experience with QGIS so I didn't found that option.

i imagine is as checking if calculated center point o district belongs to province boundaries, but maybe this is simpler task?


texts: I have a large amount of data (thousands of polygons) that I want to share with people through my site. I want them to be able to view the polygons, and due to the large amount of data, I am not worried about them stealing this data through the displayed google map. However, since the kml file must be accessible to google, a proficient coder could discover the url of the kml files and simply download the entire file. Is there any way to restrict access to these kml files only to google so no outside people can get to them. I have tried .htaccess, but I can't figure out how to restrict access to everyone but google. Any ideas?


texts:
How do I save a Standalone Table to an .mmpk file ?

I try save Standalone table to mmpk file, but i didn`t get it in mobile app

I tried get tables from this object in android app, but "GetTables" property return   zero size.

mapPackage = new MobileMapPackage(generateFileMapFoder() + "/mobile_map.mmpk");
                // load the mobile map package asynchronously
                mapPackage.loadAsync();

mapPackage.addDoneLoadingListener(new Runnable() {
                    @Override
                    public void run() {mapPackage.getMaps().get(0).getTables()})


texts: Using interactive inputs to Create Fishnet using ModelBuilder?

texts:
I'm creating a fully automatic tool in ModelBuilder, and have a problem with the Fishnet function.

The tool should work with different files, and the user only needs to choose a specific txt file (and choose the workspace) to start the process.

For this reason, I need the tool to be be able to create a fishnet extracting all the settings autonomously.
At the moment, the fishnet can extract: The Template Extent, the Cell size Width, the Cell size Height and the output feature class.

The problems are the Fishnet Origin Coordinates and the Y-axis Coordinate. In the common use of the Fishnet instrument, you just need to put the Template Extent and all these fields compiles automatically, but that doesn't happen here.


texts:
I've been following the steps in the book "An intro to using gis in marine biology, supplementary workbook seven".

I'm trying to create a presence-absence map by turning the data layer of dolphin sightings in the North Sea to a raster layer, however I'm not able to select the "PRESENCE" attribute field from the list of existing attributes in the attribute table.

Is there something I've missed?


texts: Exporting shapefile to Google Earth format (*.kmz) using QGIS?

texts: Is there a way to export a shapefile to Google Earth format (*.kmz) using QGIS?


texts: How to set map scale using Arcgis silverlight Api?

texts:
I want to set map scale manualy using ArcGIS Silverlight Api . I use silverlight api 3 and arcgis server 10 . I have 4 scale levels and I want to show the map in these scales. How can i  set these scale using Silverlight Api ?

scale :
1:288895.27714400
1:144447.63857200
1:72223.81986000
1:36111.90964300

thanks for advance


texts:
I have been tasked to make a map or some sort of analysis that analyzes the spatial relationship of a certain ethnic that has high unemployment to areas that have potential employment opportunities. Data that I have are: Current census data, employer point data, and unemployment data.

Is there any other data I might need and what arcmap extension is out there that I can use that will easily analyze the data or do I have to get ESRI Business Analyst?


texts: With QGIS, how can I create a layer from a .csv file?

texts:
I have a .csv file with ~16k points with lat and long to plot onto a map. I'm completely new to both QGIS and GIS in general and have got to a point where I have downloaded the OpenLayers and Georeferencer plugins and have a raster satellite image, but my data points are all collected off the coast of Africa at a much smaller scale than the map itself.

However, when I try to use the Georeferencer plugin it only lets me add a raster image from my local files (unless there is an option I haven't seen), yet the image I am using is from a plugin; I'm not sure if there is any way of finding it from my home directory to load it into the Georeferencer.

How do I georeference a map loaded with OpenLayers?


texts: Cannot load maps using OpenLayers plugin in QGIS

texts:
I'm currently using QGIS 2.12.1-Lyon. I have installed the OpenLayers Plugin (and updated to the latest version, v1.3.6). However, when I go to load up base maps in QGIS using this plugin (go to Web | OpenLayers plugin, then select a map) QGIS always throws a Python error.

The same error is always thrown. Some maps load, but most do not. Not at all sure why. I am still able to use the Google Hybrid base map for example, but not Google Streets.

Error dump is shown below (have changed user names and private folders for privacy).

An error has occured while executing Python code:

AttributeError: 'OpenlayersLayer' object has no attribute 'geometryType'
Traceback (most recent call last):
  File "C:/Users/matlab/.qgis2/python/plugins\BLNExporter\BLNExp.py", line 208, in EnablePlugin
    if (layer.geometryType() == QGis.Polygon) or (layer.geometryType() == QGis.Line) or (layer.geometryType() == QGis.Point):
AttributeError: 'OpenlayersLayer' object has no attribute 'geometryType'

Python version: 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]
QGIS version: 2.12.1-Lyon Lyon, 16760fd

Python Path:
C:/Users/joebloggs/.qgis2/python/plugins\processing
C:/PROGRA~1/QGISLY~1/apps/qgis/./python
C:/Users/joebloggs/.qgis2/python
C:/Users/joebloggs/.qgis2/python/plugins
C:/PROGRA~1/QGISLY~1/apps/qgis/./python/plugins
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\matplotlib-1.3.1-py2.7-win-amd64.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\requests-2.0.1-py2.7.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\python_dateutil-2.2-py2.7.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\six-1.9.0-py2.7.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\lxml-3.3.5-py2.7-win32.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\gsimporter-0.2-py2.7.egg
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs\httplib2-0.9-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\nose-1.3.3-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\tornado-4.0.1-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\backports.ssl_match_hostname-3.4.0.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\certifi-14.05.14-py2.7.egg
C:\PROGRA~1\QGISLY~1\bin\python27.zip
C:\PROGRA~1\QGISLY~1\apps\Python27\DLLs
C:\PROGRA~1\QGISLY~1\apps\Python27\lib
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\plat-win
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\lib-tk
C:\PROGRA~1\QGISLY~1\bin
C:\PROGRA~1\QGISLY~1\apps\Python27
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\PIL
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\jinja2-2.7.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\markupsafe-0.23-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\pytz-2012j-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\win32
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\win32\lib
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\Pythonwin
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\Shapely-1.2.18-py2.7-win-amd64.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\six-1.3.0-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\wx-2.8-msw-unicode
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\xlrd-0.9.2-py2.7.egg
C:\PROGRA~1\QGISLY~1\apps\Python27\lib\site-packages\xlwt-0.7.5-py2.7.egg
C:/Users/joebloggs/.qgis2//python
C:\Users\joebloggs\.qgis2\python\plugins\opengeo\ext-libs
C:\PROGRA~1\QGISLY~1\apps\qgis\python\plugins\fTools\tools
C:/Users/joebloggs/Documents/RandomFolder/RandomSubFolder


texts: Preserving styles when going from KML to WMS

texts:
I'm currently losing styling information somewhere along my chain of applications:

KML -> ArcGIS Layer (via ArcGIS Desktop)
ArcGIS .gdb -> PostGIS via ogr2ogr
PostGIS -> GeoServer
GeoServer -> Leaflet via WMS service

How might I preserve styling information when serving data originally stored as a KML via WMS? I have found no way to export styles from a KML into SLDs.

If this can natively be done through a different type of GIS Server, please advise.


texts: How to solve a filename conflict problem?

texts:
I use a procedure to remove and delete a feature layer with a certain name. Then, using another procedure which is called during the run time of the program (as the remove and delete procedure as well) I try to create a new feature layer which needs to have the same filename as that removed-deleted before. Unfortunatley, the program stops in the code line that tries to create the new feature layer. It is definetely because of the conflict between the previous and the new filename which are the same. How can I solve this problem?
Any help please?


texts: How to manage the CRS of OpenLayers plugin in QGIS?

texts:
I am working with OpenLayers plugin in QGIS 1.8 but I am not sure on how should I manage the CRS.

The default CRS for Google/Bing layers is WGS84/Pseudo Mercator (EPGS:3857) but I am not sure on the CRS of the new layer I want to work on.


texts: Creating hotlinks only on objects that have value using MapBasic?

texts:
In mapbasic I am creating hotlinks with the code below. However when there is no value in the column it creates a hyperlink that doesn't work and this is frustrating for users. Is it possible to only have the hyperlink on records with a value? or at least for those without a value the alias should say "No AsCon found"

Map From temptable
mapWinID = FrontWindow()
Set Map Window mapWinID Layer 1 Activate Using "https://server.com/docs/~D"+"As_con_InfoExpert_ID" On Labels Objects Relative Path Off Enable On
Metadata Table temptable SetKey "\ActiveObject\Expr" To ("""https://server.com/docs/~D""+As_con_InfoExpert_ID")
Metadata Table temptable SetKey "\ActiveObject\Mode" To "HOTLINK_MODE_BOTH"
Metadata Table temptable SetKey "\ActiveObject\Relative" To "OFF"
Metadata Table temptable SetKey "\ActiveObject\Enabled" To "TRUE"
Metadata Table temptable SetKey "\ActiveObject\Alias" To ("Open in Infoexpert")


texts: I tried to find QGIS plugin "Load raster to PostGIS" but I could not find anything in Plugins section in QGIS. Is there any external resource where is possible download this plugin and after that install in QGIS 3.4?


texts:
When I start the GPS (in USB), I look at the data comming and after a while, first reasonable data is $GPGSV sentence (see http://aprs.gids.nl/nmea/#gsa), that contains azimuth and elevation of satellites. Only after some time, real coordinates start to come in another sentence.

I thought first, it may be some world coordinates in $GPGSV sentence, but the values are in the ranges expected for azimuths (0-360) and elevations (0-90).

How could my GPS know satellite's elevations (and azimuths) before knowing its own actual location?

Do I miss something? Is the satellite signal similar to a VOR signal?


texts:
I've been playing with Bing Maps + Polymaps recently.

It has become apparent that Bing is the only mapping resource available to Polymaps that includes satellite imagery. Is this true?

I am trying to accomplish drawing many polygons of different sizes on a map, and filling them with different colors depending on what they represent. Furthermore, I want each of these polygons to be clickable.

I understand that to load geographical data onto a map with Polymaps, you either use .features() or .url() as part of the geoJSON object.

I have .features() working, but by using this function the polygon data must be embedded into the javascript file itself - as in:

map.add(po.geoJson()
        .features([
    {
      "geometry": {
        "coordinates": [[
[[-63.599925935268395,44.64018531837119],[-63.599789142608664,44.640219671024646],[-63.59973818063736,44.6401452402498],[-63.599877655506134,44.64010707058464],[-63.599925935268395,44.64018150140844]],

[[-63.59986692667008,44.639874235083376],[-63.59975695610047,44.63991431339196],[-63.599697947502136,44.63981507181532],[-63.5998159646988,44.639782627416906],[-63.59986424446105,44.6398723265918]]
        ]],
        "type": "MultiPolygon"
      },
      "type": "Feature",
      "id": "beat:05X",
      "properties": {"name": "05X"}
    }

        ]));

I want this data to actually come from a database and not lie statically within the javascript itself - which tells me I should use the .url() method instead - correct? And if I use the .url() method, I must pass a valid URL template that returns geoJSON... How do I go about this with my own private data? I'm familiar with PHP and mySQL - must I use a spatial database system instead?

And how does a map server like MapServer, GeoServer, and others fit into this? Could I not set up a PHP script to read these templates from the $_GET variable, query a database and return the proper geoJSON string?

I really just need some inspiration, a few pointers to set me on my way and kill some of the doubts I have. I've found myself just digging around the stackexchange and googling everything from map servers, to geo spatial databases - you name it, I've probably read about it, but I'm struggling putting the pieces together.

Long winded question but once these things are clarified, I should be in a much better position.

We can remove polymaps from the equation if there are substitutes.

In summary, I'm looking to take a large number of polygons from a database (PostGIS?) of all different shapes, sizes, colors, and transparencies, and plot them on a map. That's essentially it! Polymaps has vector graphics + I thought getting this on the go would be trivial given the above code... we can scrap it if you have other ideas. What is the best solution? My only business requirement is that Satellite Imagery be enabled...Which means we're working with Bing or Google.


texts: Calculate aspect for stream segments in R from 2D linestring and DEM

texts:
I would like to calculate the general aspect for stream segments in R. The stream network has been derived from the DEM and is in 2D. I have split the stream segments into equal parts of 100m (or less) and I have the DEM (1m resolution). While the aspect may change (slightly) along each segment, I am interested in the overall or dominant aspect for that segment based on the hill slope from start to finish. [I can see it might possible to use the DEM to create a 3D linestring but I don't know how]
When I use a basic tool to calculate aspect for each cell/pixel of a DEM raster, it shows opposite facing-slopes on either side of the stream (logically). Here, I want the aspect of the stream, i.e. of the general slope/aspect of the hill along which the stream flows. If the stream is flowing along a hill that is north-facing, the aspect that I am hoping to obtain is North (or 0). I don't want the direction of flow of the stream (from x1,x2 to y1,y2) although it may overlap given that the stream will probably flow in the most direct route downhill.
I am aware of tools in QGIS and ArcGIS but I would like to make as many of my steps within R as possible (Calculate slope of line segments with QGIS ; Workflow for determining stream gradient?)
I had initially asked a question about slope and aspect but only got an answer for slope so I'm asking it specifically about aspect here.


texts: Postgis/Postgresql - Interpolate equadistant points based on weighted attribution similar to QGIS Interpolation(Heatmap Kernal Density Estimation)

texts:
I have equadistant vector points, that are all 400m apart but only cover certain areas of a map that have a value of (1-5). In QGIS I can run the Kernel Density Estimation with my vector points to select a distance, pixel rows, and then weight based on the (1-5) values, 1 best/5 worst to get an interpolated distribution. Here is the before and after. What I'm attempting to do is write a script that can essentially replicate what QGIS is doing but using POSTGIS functions to have an input of vector points and outputs of vector polygons such as this. Using the Polygonize Raster to Vector after Kernel Density is run and then put the new polygon shapes into a Graduated DN into 5 bins of equal interval:

Input:

Kernal Density Params & Output:

Then toward the end of the process I will then dissolve the similar bins to have a seemless display such as this:


texts: How to detect green areas in Google Earth imagery using R?

texts:
I want to detect urban green areas in Google Earth imagery, using R. I've already downloaded contiguous images and assembled them into a larger image.

The above PNG image has 800x769 pixels, but the original has 6400x6150. When I open it in R, using png library (and readPNG function), I get a matrix with three layers, containing the red, green and blue channels (same as with Gimp).

I want to know if it's feasible to extract only the green areas (areas with native vegetation) - so that later I can transform them into shapefile polygons or raster masks. Note that some areas with native vegetation are more brown than green (savanna type vegetation), and a big green area on the northern half of the map correspond to a lake.

I know that's usually done with non-visible bands of satellite spectrum, but I want to know if I can do this with vanilla Google Earth.


texts: Qgis 1.9 commitChanges() and Redraw a temporary layer

texts:
Using QGIS 1.9 and writing a python plugin.  I have a function that connects points into polylines and draws them as a temporary, in memory layer.  Membership in a polyline is determined by a field in the attribute table (postgis enabled table) of the point.

The tool programatically sets up an edit session and alters the correct column in the point attribute table to alter the polylineID.

pseudocode:

layer.startEditing()

...Make edits to the features (not a data provider)...

layer.commitChanges()
remove in-memory layer
redraw in-memory layer by iterating over the point table to building polylines

The redraw fails to see the edits that have occurred.  This is after calling layer.commitChanges().  What is the proper method to make an edit and propagate it so that an immediately derived layer sees the changes?

Tested with data in geojson on disk and via a view on a postgis table.


texts: If statement in ModelBuilder to check if field exists

texts:
I'm trying to create the following in ModelBuilder:

Set input feature class
  if field1 exists
    execute the processes in branch1
  else
    add field1
    execute the processes in branch2
    merge branch

While I have no problem with the processes in the two branches, the IF check at the start doesn't seem to be working the way I want it to work. I am using this code in the Calculate Value code block

def ifFieldExists(Input_FC, SUMMARY):
  import arcpy
  fields = arcpy.ListFields(Input_FC)
  for f in fields:
    if f ==" SUMMARY":
      return "true"
    else:
      return "false"

and expression

ifFieldExists(r"%Input_FC%", "%SUMMARY%")

This is the precondition for branch1, so if it returns true, the field exists so branch1 executes. The other branch's precondition returns the reverse. When I run it on an input that does contain the field, the precondition returns false instead of true, so the wrong branch executes.


texts: Selecting attributes that contain "%"-symbol in Field Calculator

texts:
I have a Text field and it contains values that have a "%"-symbol and "°"-symbol. Example:
1.2% 3.2° 2.3% -1.2°

I don't know how to select just the values that have the "%"-symbol because Field Calculator thinks its a Wild Card. Any ideas?


texts: Updating raster value to reflect total number of points within cell using Python?

texts:
I have to do a task in Python/ArcPy/GDAL/whatever because I'm trying to automate this process for a hundred or so different raster/shapefile combos.

This is what I need to do:

I have an empty raster and a shapefile with thousands of points. I would like to update the value of each cell/pixel in the raster to reflect the number of times a unique point lands within that cell.

At first I was thinking something like this:

for cell in raster:
    for point in shapefile:
         if point intersects with (/is contained in) cell: cell_value + = 1

However, I don't know how to implement this because I don't see how to intersect a raster and a polygon. Converting the raster into polygons is out of the question because the rasters are very large with a lot of cells.

I'm not married to the idea of looping through both the individual cells and points because I feel like it's unnecessary, so if there is an ArcGIS tool(s) that I could manipulate, and/or another method in python I could use to reach my goal, I'm more than happy to try that out.


texts:
I have an attribute table of xy points that I have plotted on a map, and a shapefile that divides the map into 46 different areas.

I want to add a field/variable to the table containing the xy points, that states which area of the map the point is located.

How do I go about doing this?


texts: Does QGIS have an autocomplete polygon option like the one in ArcGIS?

texts: I am trying to do digitizing in QGIS, but I do not see an option that I have seen in ArcGIS - autocomplete polygon. I would like to make two adjacent polygons topologically correct, and can not find such an option in QGIS.


texts: What language is used in QGIS 3 label expressions?

texts:
What language is used in QGIS 3.x's label expressions?

Initially, I thought it's python. But when I tried to concatenate a string and a number field as follows (in Layer property Labels -> Single Labels):

name + str(id)

I got an error about Function is not known (I guess str()).

Is the language here python, PostSQL or something else?


texts: Does anyone know if there is any open source data downloads for Northern Ireland, seems not as straight forward as rest of UK?


texts: Certain MultiPolygons cause “BOOM! Could not generate outside point!” in PostGIS?

texts:
I've just posted this at StackOverflow and someone suggested I try here.
Original: https://stackoverflow.com/q/5863859/246265

I'm trying to represent a rectangular area which crosses 180 degrees longitude. For more background see https://stackoverflow.com/q/5737506/246265

Here's my test case:

from django.contrib.gis.geos import Polygon, MultiPolygon
from my_project.my_app.models import Photo

a = Polygon.from_bbox((30, -80, 180, 80))  # the part to the east of 180
b = Polygon.from_bbox((-180, -80, -160, 80))  # a part to the west of 180
c = Polygon.from_bbox((-180, -80, -100, 80))  # a larger part to the west of 180

ok   = MultiPolygon(a,b)
ok2  = MultiPolygon(c)
boom = MultiPolygon(a,c)

# This works
Photo.objects.filter(location__coveredby=ok)[:1]
# This also works so c is ok
Photo.objects.filter(location__coveredby=ok2)[:1]
# This gives "BOOM! Could not generate outside point!"
Photo.objects.filter(location__coveredby=boom)[:1]

# splitting c doesn't help
c1 = Polygon.from_bbox((-180, -80, -140, 80))
c2 = Polygon.from_bbox((-140, -80, -100, 80))
test = MultiPolygon(a,c1,c2)
Photo.objects.filter(location__coveredby=test)[:1]
# BOOM! Could not generate outside point!

By changing the numbers I can make this error come and go.
(-180, -80, x, 80) works where x <= -140 for example. For every number there is a threshold like this but I can't find a pattern. For boxes with the same area, some will work and others not. For boxes with the same width some will work and others not.

I can look at the SQL being generated but the areas are represented in binary (EWKB) and I'm not sure how to read it.

Can anyone explain this?


texts: Getting error when trying to use IPyConsole Plugin in QGIS 2.18.16

texts:
I am trying to use the IPyConsole Plugin and the error I receive tells me I have to install IPython 3.1.0 or Jupyter 1.0.0. So I followed the link they gave me and installed Anaconda 2.17 and then 3 but still receive the same error. I am using Windows 7 Enterprise.
The exception error that is popping up is:

NO module name IPYthon.qt.console.rich_ipython_widget.


texts: How do I add a text file to ArcMap using ArcPy?

texts:
I'm specifically referring to the output of the Compare Table (Data Management) tool, which is a comma delimited text file.

When you add the text file manually, you can actually view it like any other table, by selecting "View Attribute Table," even though it's a just text file.

Previously, I've written code that adds the output of other tools to the map document, but these have always been feature class outputs, not text file outputs. Here's how I've added feature class outputs to an mxd in the past:

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
fc = pathToMyFeatureClass
arcpy.MakeFeatureLayer_management(fc, "fc_lyr")
addLayer = arcpy.mapping.Layer("fc_lyr")
arcpy.mapping.AddLayer(df, addLayer)

So I'm looking for something similar for the text file, and I'm not finding anything. Seems simple enough. I thought maybe something like arcpy.MakeTableView_management() might do the trick, but when I added that to my code, I got an error that my parameters were invalid, that the output of my Compare Table tool did not exist or is not supported, and since I know the output exists, then it must not be supported.

Per below comment, here are some screenshots. The attributes are sort of nonsensical to protect privacy. First screenshot is of the output text file, and the second is what it looks like once you view it in ArcMap.


texts: What is the precision of SELECT DISTINCT on PostGIS geometry column?

texts:
I wonder what the precision of the SELECT DISTINCT operator is on a PostGIS geometry. On my system, the following query gives me a count of 5, which means that the inserted points are considered equal if they differ by less than 1e-5 and I am not sure if that is a feature of PostGIS, a problem of my installation or a bug.

Does anyone know if that is the expected behavior?

CREATE TEMP TABLE test (geom geometry);
INSERT INTO test
    VALUES
        (St_GeomFromText('POINT (0.1 0.1)')),
        (St_GeomFromText('POINT (0.001 0.001)')),
        (St_GeomFromText('POINT (0.0001 0.0001)')),
        (St_GeomFromText('POINT (0.00001 0.00001)')),
        (St_GeomFromText('POINT (0.000001 0.000001)')),
        (St_GeomFromText('POINT (0.0000001 0.0000001)')),
        (St_GeomFromText('POINT (0.00000001 0.00000001)')),
        (St_GeomFromText('POINT (0.000000001 0.000000001)'));

SELECT COUNT(*) FROM (SELECT DISTINCT geom FROM test) AS test;

 count
-------
     5
(1 row)

I am using:

$ psql --version
psql (PostgreSQL) 9.3.1

and

SELECT PostGIS_full_version();
----------------------------------------------------------------------------------------------------------------------------------------------------------------------
POSTGIS="2.1.1 r12113" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.10.1, released 2013/08/26" LIBXML="2.7.3" LIBJSON="UNKNOWN" RASTER

on OSX 10.9


texts: Removing overlapping lines from two line features

texts:
I have a street network and another layer with a network of pipes, both as lines features. Some of the pipes overlap the street network.

I only want those pipes that follow along the streets. Means those which are equal to the streets or in other words remove all pipes which overlap the street network.

How can I do this with QGIS?


texts: Populate_Geometry_Columns does not work with materialized view

texts:
I have created materialized view with Points as ... ST_SetSRID(ST_MakePoint( longitude,latitude), 4326) as geom... column data. The table kind of works, I can create spatial index which works etc, but the geometry_columns view shows srid=0 and type=GEOMETRY. This is major issue for QGIS which seems to try to full scan the table to find out true geometry type. And that takes too much time as my view is quite large (over 500M rows).

Next I tried to fix the metadata with select Populate_Geometry_Columns('public.mymaterializedview'::regclass) , but this returns immediately '0' instead of '1' what I normally get with real tables. And the geometry_columns does not get changed.

I tried also setting column type with ALTER TABLE public.mymaterializedview ALTER COLUMN geom type geometry(Point, 4326) USING ST_SetSRID(geom, 4326); but this gives [42809] ERROR: "mymaterializedview" is not a table, composite type, or foreign table which is not very helpful for me.

The DB is GCP managed Cloud SQL DB, i.e. Postgres 9.6 and PostGIS full version is currently POSTGIS="2.3.0 r15146" GEOS="3.5.0-CAPI-1.9.0 4392" PROJ="Rel. 6.0.0, March 1st, 2019" GDAL="GDAL 2.4.0dev, released 2018/99/99" LIBXML="2.9.1" LIBJSON="0.12.1" RASTER


texts: How can I get a KML/KMZ file to display on Google Maps without a public facing web server?


texts: OGR: How to save layer from AttributeFilter to a shapefile?

texts:
I select a subset of features from my shapefile like so:

node_lyr = nodes.GetLayer()
q_str = "%s = %i" % (field, value)
node_lyr.SetAttributeFilter(q_str)

What is the best way to make a shapefile from this subset?


texts: How to understand the which GIS server is used?

texts:
I want to understand how I can learn the GIS server type of a certain link. In instance I want to know which type of server is being used for http://barajlar.dsi.gov.tr/ web site.

Related server links are http://tilesrv01.dsi.gov.tr/ and spatialsrv01.dsi.gov.tr. Are they using GeoServer, ArcGIS, Mapnik, ..vs?


texts: How to calculate the percentage of a Line through Polygons?

texts:
I have a straight line which runs through a bunch of polygons, which don't have spacing between them (border to border).

Is there an easy way with QGIS to calculate which polygons it passes and by what percentage of the total distance its within these polygons?

I'm using Ubuntu, QGIS 1.8.


texts: Error executing tool.: Layer "layer1" Unable to connect to map server

texts:
I am using the Configurable Javascript Viewer, and have only changed a couple of the configuration settings.

One change I made was to use my own Export Web Map Task.

In viewer.js, I changed:

printTaskURL: 'http://maps.evansvillegis.com/arcgis_server/rest/services/Utilities/PrintingTools/GPServer/Export%20Web%20Map%20Task'

I then saw this note in the readme.md in the proxy directory:

  You will more than likely need a proxy page for printing and other large cross domain requests.

So, I set that up using the provided proxypage_net.zip files, proxy.ashx and proxy.config.

In proxy.config, I added:

<serverUrl url="http://maps.evansvillegis.com"
           matchAll="true"></serverUrl>

Then I made sure my viewer.js file pointed to the current proxyUrl:

esriConfig.defaults.io.proxyUrl = 'proxypage_net/proxy.ashx';

However, when I try to print, I am getting the following error in the response:

{"error":{"code":400,"message":"Unable to complete operation.","details":["Error executing tool.: Layer \"layer1\": Unable to connect to map server at http://maps.evansvillegis.com/arcgis_server/rest/services/BASEMAPS/BASEMAP_GOOGLE_WGS84_NO_FOOTPRINTS_BW/MapServer.\nFailed to execute (Export Web Map).\nFailed to execute (Export Web Map Task)."]}}

I have seen related posts:

Secure Print Service Failed to Execute

ArcGIS Online printing Unable to connect to map server?

However, the provided information (mind you, limited) was not much help.

The 400 error - is this because my proxy is not setup correctly, or is there something wrong with my PrintTask service?

This is the default PrintTask that is configured with ArcGIS Server 10.2.1.  The only other change I have made, is ensuring IIS users and the ArcGIS Server Account can access the ArcGIS Server directories.

Viewer: http://maps.evansvillegis.com/viewer


texts: QGIS Python API set label color based on attribute field

texts:
I'm trying to set the text label color and size based on the attribute table. I got the size setting correctly but not the color.
My current result looks like this (The symbol has been set to be transparent), and I want to give different colors to different text labels.

The attribute table looks like this:

The code is as following:
vlayer = QgsProject.instance().mapLayersByName(layer_name)[0]
vlayer.renderer().symbol().setOpacity(0.0)

pc=QgsPropertyCollection('Text')
color_prop=QgsProperty()
color_prop.setField("color")

size_prop = QgsProperty()
size_prop.setField("size")
pc.setProperty(QgsPalLayerSettings.Color,color_prop)
pc.setProperty(QgsPalLayerSettings.Size,size_prop)

label_settings = QgsPalLayerSettings()
text_format = QgsTextFormat()

label_settings.setDataDefinedProperties(pc)

label_settings.fieldName = 'name'
label_settings.setFormat(text_format)

vlayer.setLabeling(QgsVectorLayerSimpleLabeling(label_settings))
vlayer.setLabelsEnabled(True)
vlayer.triggerRepaint()

I wonder where is wrong with the color setting in the code? My best guess is that the color field in the attribute table is not what QGIS requires, but I'm not sure what should be the correct format.


texts: Earth Engine: Flatten (?) list of filters for ee.Filter.or

texts:
Is there a more elegant way than calling each filter in the list? It looks like ee.Filter.or wants separate arguments, not a list.

Using var NEfilter = ee.Filter.or(filterlist) gives the error:

Filter.or, argument 'filters': Invalid type. Expected: List< Filter >. Actual: List< List< Filter >>

https://code.earthengine.google.com/dffaf31424c0b5cbf2fd69187523eae0

var AllStates    = ee.FeatureCollection('ft:1fRY18cjsHzDgGiJiS2nnpUU3v9JPDc2HNaR7Xk8');   // U.S. States
print (AllStates);

var NameFcn = function(name) {
  return ee.Filter.eq('Name', name);
  };

var NE_list = ["Massachusetts", "Connecticut", "Rhode Island", "Maine", "Vermont", "New Hampshire"];
var filterlist = NE_list.map(NameFcn)

print(filterlist)

var NEfilter = ee.Filter.or(filterlist[0],filterlist[1],filterlist[2],filterlist[3],filterlist[4],filterlist[5]);

var NEStates = AllStates.filter(NEfilter);
print(NEStates)


texts: List ArcGIS geodatabase feature classes on Python console using ListFeatureClass function

texts:
I am trying to create a new script with a hard-coded geodatabase filename, that lists the feature classes it contains to the console. I know there are several ways to do this, but I decided to to use ListFeatureClasses function.
 I barely started studying Python, so I know very basic stuff, but I would like to get help with this.

I tried the following code, but it does not print anything and does not even give me any error message. I opened up ArcGIS Pro, added the geodatabase to my map, and wrote this code in the built-in Python console. Is there any step I am missing or something wrong with my code?

import arcpy
import os

arcpy.env.workspace = "X:\311\Obtaining GIS Data\TaxParcels.gdb"

datasets = arcpy.ListDatasets(feature_type='feature')
datasets = [''] + datasets if datasets is not None else []

for ds in datasets:
    for fc in arcpy.ListFeatureClasses(feature_dataset=ds):
        path = os.path.join(arcpy.env.workspace, ds, fc)
        print(fc)


texts:
How do I create oriented centerlines in QGIS as in the following image:

Input data is a shapefile of irregular polygons as in the image below and available on the following link:https://drive.google.com/drive/folders/1fAS3BDDis6z154L0nJVelWJoKzU04lFN?usp=sharing

Figure 1: Shapefile of irregular polygons
Secondly, Minimum bounding box tool was used to create a bounding box ('perfect rectangles') around each polygon as in the image below and available on the following link:https://drive.google.com/drive/folders/1Tr3umzfsE2dBVna9J1tYE89eoe3Vc0gR?usp=sharing.
QGIS > Processing Toolbox > Vector Geometry > Minimum bounding box

I need to create oriented centerlines to produce rows of lines as in the example image below:


texts:
What I am trying to do is create a table contains intersections at the location where each road intersects Roadname 1 & Roadname 2 as well as the geometry spatial data for the intersection ST_Intersection(a.geom, b.geom)

But it  seems it doesn't work out l. I am pretty new to PostgreSQL, and I am wondering if there is a way to do so? The following is what I have so far.

Create Table roads(
    intid_ integer,
    muni text,
    alias_ text,
    roadnames text,
    geom geometry);

Insert into roads (roadnames, geom, intid_, alias_)
SELECT
    trim(student.sc_e911slrn.prefix) || trim(student.sc_e911slrn.name) || trim(student.sc_e911slrn.suffix),
    geom, objectid, alias
FROM student.sc_e911slrn;

--------------------------------------------------------------------

CREATE type student.inter AS (
     roadnames text,
        intid_ int,
        alias_ text,
        muni text,
        ST_Intersection geometry,
        geom geometry);

--< Cursor >--------------------------------------------------------
Create function create_intersections()
    /*returns table student.intersect_roads As $method1$*/
    RETURNS table (
        roadnames text,
        intid_ int,
        alias_ text,
        muni text,
        ST_Intersection geometry,
        geom geometry) As $$

Declare
    return_intersect inter;
    rec_road Record;
    cur_road Cursor()
        For Select intid_, geom, roadnames
        From student.roads;

Begin
    Open cur_road();

    Loop
        fetch cur_road into rec_road;
        exit when not found;

        If ST_Touches(rec_road.geom, student.roads.geom) AND
            rec_road.geom != student.roads.geom AND
            rec_road.intid_ != student.roads.intid_

        Then Update student.intersect_road
             Set st_intersection = ST_Intersection (rec_road.geom, student.roads.geom),
                 intid_ = rec_road.intid_,
                 intersections = rec_road.roadnames|| ' & ' || student.raods.roadnames,
                 alias_ = rec_road.alias_,
                 /*muni = */
                 polyline = ST_Union(rec_road.geom, student.roads.geom)
                 RETURNING * into return_intersect;
                 RETURN NEXT;
        End if;
    End Loop;
    Close cur_road;

    RETURN;

End;
$$ Language plpgsql;

--< return the table>----------------------------------------------
SELECT create_intersections();

The Roadname 1 & Roadname 2 should be generated from roadnames column from ROADS table, and the geometry for the intersection ST_Intersection(a.geom, b.geom) should be generated from the geom column from ROADS table.
When I run it, it returns no result.

The table ROADS is something like this road table,(geometry view):

Intid_     Muni    Alias_           Roadnames            geom
24581722   NULL    CR 52       Coopers Falls road        ...
24581723   NULL    CR 44       Rama road                 ...

Expected output should be similar to this one, but we need other fields to be filled as well.--->table (geometry view)

I also tried

CREATE FUNCTION create_intersections()
     UPDATE student.intersect_roads AS $$
Declare
Begin
End
$$

Which is not working as well.


texts: NameError: global name is not defined when performing task with FlowMapper in QGIS

texts:
I'm using QGIS 2.8 and I'm trying to apply the FlowMapper plugin to create a visualization. I aim to follow the instructions in this blog post and this GIS.SE question, but something seem be going wrong.
The error is as follows:

2020-11-04T19:45:59   1   Traceback (most recent call last):
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/FlowMapper/flowmapper.py",
line 2655, in run
SaveDirectory = SaveShpName             NameError: global name 'SaveShpName' is not defined

I'm not sure what to do about this. So far, I've tried changing the input files by either switching to spaces instead of tabs between the data columns. That did not help though.
Added: Here are the files I've used. First, there is the locations of the cities:

Then there's the labels associated with them:

And finally, the incidence matrix:

All columns are separated by spaces. The paths to the respective files are:

/Users/maxmuller100/Documents/Data_Project/stud-steden-locaties.txt

and

/Users/maxmuller100/Documents/Data_Project/stud-steden-labels.txt

and

/Users/maxmuller100/Documents/Data_Project/stud-steden-matrix.txt


texts: Setting output parameter in Python script tool to be PDF?

texts:
I need to allow a user to input a path and name where they would like to save their new map book PDF on their computer/drive.

I have created a script where the user can input their mxd with data driven pages, input their title map, append the two PDF documents and save the appended PDF, but it's all hard coded and I'd like to create a script where they can select where to save their new appended PDF when they open it in the toolbox.

I need the script to be dynamic. I'm stumbling on setting up the parameters. I'm unsure what the Output PDF data type should be set as in the toolbox. None of the selections seem to fit.

Aside from "Output" in the "Parameter Properties" should I change anything else?


texts: Error sending request to OpenRouteService API from a set of client addresses to a foodbank site

texts:
I am sending a request to OpenRouteService API from a set of client addresses to a foodbank site:

I am passing coordinates in the form ((x1,x2),(y1,y2)) where x coordinates are from a large csv file, and the y coordinate is the extra appended coordinate.

I sent a successful request (for a single coordinate pair) to openrouteservice in the form

cords = ((8.34234,48.23424),(8.34423,48.26424))

I am now trying to send a request for a large batch of coordinate pairs.

import csv
import openrouteservice
client = openrouteservice.Client(key='myAPIkey')

with open ('ClientCoordinates.csv' , 'r') as csv_coords:
   csv_reader = csv.reader(csv_coords)
   EXTRA_COORD = [53.464774, 2.23179]
   csv_coords = []
   for line in csv_reader:
       csv_coords.append([line, EXTRA_COORD])

coords = (csv_coords)
routes = client.directions(coords)

print(routes)

This returns ApiError:

  400 ({'error': {'code': 2002, 'message': "Parameter 'coordinates' has
  incorrect value or format."}, 'info': {'engine': {'version': '5.0.1',
  'build_date': '2019-05-29T14:58:37Z'}, 'timestamp': 1567170534336}})

Perhaps my use of brackets is incorrect, as OpenRouteService requires coordinate inputs in ().


texts: Temporal coverages dynamically reloaded from data store in GeoServer

texts:
I'm pretty new to GeoServer. I would like to serve results of a simulation that concerns geospatial data as a Web Coverage Service (it has to be OGC compliant). The individual steps of the simulation (similar to movie frames) should be uploaded to the server. Currently the results are simply series of files with some basic metadata. Client may want to ask for a specific period of time that should be returned (specific set of frames).

I would like to write my own data store format (GeoServer extension) for reading simulation results that would be kept in some directory. However, I see 3 problems with doing this with the GeoServer:

There is only a 1-1 relation between data store and layer. That is I would have to create separate data stores for each simulation. Because each simulation is of the same kind it would be more reasonable to keep them in one store.
The layers do not refresh themselves at runtime. That is, if a new data appears in a data store then it is not visible by the layer. The layer needs to be recreated.
I am worried about the support of time-querying of coverages or multidimensional coverages. I read somewhere that those are not supported well by the GeoServer.

Another approach is to create layers and data stores dynamically using REST interface when new simulation results are available. This doesn't look like a proper way, though.

Am I right about above issues? Is there any convenient way for doing this with GeoServer?

I am using GeoServer 2.2-beta.

EDIT:

We eventually ended up implementing our own simple WCS service that meets our needs. However, because the ImageMosaic plugin did work when it comes to series of time-related data (the problem was we still had to create new layers manually) I accept Andrea Aime's answer.


texts: We are using MapInfo in our organisation to maintain corporate GIS datasets (in TAB files). We've come to realise the importance of maintaining a metadata file for our GIS layers. Are there any free/open-source metadata managers out there that can help us achieve the desired objectives?


texts: Application of a Piecewise Logistic Function to Extract Seasonality Metrics

texts:
I am working with 16-day MODIS EVI (satellite) data and I want to fit a Piecewise Logistic Function through my 23 EVI data values. The following formula is for the Piecewise Logistic Function:

y(t)= [c / (1 + e^(a+bt))] + EVIb

where:
t is the day of year
a is vegetation growth time
b is rate of leaf development
c is the amplitude
bEVI is the background EVI value [constant]

For some reason there is very little clarification about each component of the equation in most papers. I am assuming that the bEVI value is a constant but am supposing that one or two other components are constants too?

I have little experience with this function, but from the literature I have read, is applied separately to the left and then subsequently to the right of this bell-shaped curve - split at the maximum y value. I am hoping someone here with a stronger mathematical background and experience using the piecewise logistic function would be able to advise me on which component of the equation above should be a constant and which should be a variable.

I tried embedding the original journal paper that I have on my OneDrive but can't yet seem to figure out how to do it. I have, however, included the link here "Monitoring Vegetation Phenology using MODIS (Zhang et al., 2003)". You may require access to Science Direct but please do let me know if you have any problems with that. The piecewise logistic regression is mentioned on page 472, section 2.

I contacted the original author regarding the equation and he clarified that the function if applied separately for the left and right side of the curve, split at the maximum point, which could most likely be somewhere in the middle of the curve.

The series of days of the year are as follows:
1
17
33
49
65
81
97
113
129
145
161
177
193
209
225
241
257
273
289
305
321
337
353

Here is an example EVI dataset I currently have for one year of 23 data points:
0.425
0.3899
0.2417
0.4724
0.4648
0.226
0.245
0.2509
0.2741
0.3711
0.3914
0.3985
0.4433
0.397
0.3329
0.3578
0.2741
0.2912
0.3053
0.2995
0.0959
0.3287
-0.3


texts: PostGIS - geocode and address standardization

texts:
Can you please help me ( I am not professional, If it is possible some easy explanation) how to proceed with address standardization and geocoding with Open Street Maps data (outside the US) in PostGIS?

I've add osm data into PostGIS over osm2pgsql, but I edit style file for columns like street, country code, city ... So I think I have all important informations about address. What should I do next? I need some steps how to proceed, and I next study required technologies.

I am working on Windows.


texts: QGIS atlas generation - batch process for many atlases

texts:
I have many different QGIS projects, each one is designed to create one particular atlas or map book in pdf format (one different book per every QGIS project).
Month by month I must actualize each atlas, so I have to open each QGIS project and generate the corresponding pdf book. Because this process take a lot of productive time I want to find a way to actualize those atlases in a quick way, so my question is:
Is there a way to create a batch so the pdf actualization process can be automated?


texts:
I am trying to migrate a dataset from one fileGDB to another using the Append tool.

When running a few features get copied over however the rest get the following error:

“WARNING 000594”: Input feature #: Unexpected operation”

Same spatial reference.

To be clear this is an issue with the append tool as it does not work even if being run manually with the tool. I have also tried playing around with the following setting with no luck:

arcpy.env.overwriteOutput = True
arcpy.env.outputMFlag = "Same As Input"
arcpy.env.outputZFlag = "Same As Input"
#arcpy.env.XYDomain = "-180 -90 180 90"
arcpy.env.extent = "MAXOF"

# SPATIAL REFERENCE
# 26918 = UTM NAD83 ZONE 18N
spatial_ref = arcpy.SpatialReference(26918)
arcpy.env.inputCoordinateSystem = spatial_ref
arcpy.env.outputCoordinateSystem = spatial_ref
#arcpy.env.geographicTransformations = "NAD_1927_To_NAD_1983_NADCON"

# XY Tolerance
arcpy.env.XYResolution = 0.000001

def get_field_mappings(fc_in, dictionary):
field_mappings = arcpy.FieldMappings()
field_mappings.addTable(fc_in)
print ("\nField Mapping: {0}\n".format(fc_in))

for input, output in dictionary.iteritems():
    print 'Field Map: "{0}" -> "{1}"'.format(input, output)
    field_map = arcpy.FieldMap()
    field_map.addInputField(fc_in, input)
    field = field_map.outputField
    field.name = output
    field_map.outputField = field
    field_mappings.addFieldMap(field_map)
    del field, field_map

return field_mappings

def map_and_append_features(fc_in, fc_out, field_dictionary):
    field_mappings = get_field_mappings(fc_in, field_dictionary)
    arcpy.management.Append(fc_in, fc_out, "NO_TEST", field_mappings)


texts:
i'm trying to load a geojson file on a leaflet map. here is the code :

body :

<!DOCTYPE html>

    <title>map</title>
    <link href="leaflet.css" type="text/css" rel="stylesheet">
    <link href="index.css" type="text/css" rel="stylesheet">

javascript :

var map = L.map('map',{
center:[5,28],
zoom: 3,
minZoom: 2,
maxZoom: 18,
});

L.tileLayer('http://{s}.tiles.mapbox.com/v3/ianmule.bg2v5cdi/{z}/{x}/{y}.png',{attribution: "Mapbox"}).addTo(map);

function basement (feature,layer){
layer.bindPopup("SURPRISE")
}

    var c2 = L.geoJson(loc_infra, {
            style:
            {
                weight: 4,
                color: "red",
                opacity: 1,
                fillColor: "red",
                fillOpacity: 0.5
            },
            onEachFeature: basement

        }).addTo(map);

and a little CSS :

html, body, #map {
height: 100%;
width:100%;
font-family:  : sans-serif;
}

and the geojson i'm trying to launch  :

var loc_infra =
{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [........

but my layer is not visible on the map, and i can't understand why !


texts: Calculating shape.length@{Unit} on Unprojected CRS

texts:
Using the CalculateField_management tool, one can specify the unit of measurement when calculating shape length, for example Shape.length@FEET.  I was worried about doing this on data in an unprojected CRS.  Related to that, Esri Help gives the following note:

  If the data is stored in a geographic coordinate system and a linear unit (for example, feet) is supplied, the length calculation will be converted using a geodesic algorithm.
  --Esri

My question is, what does this note mean and what things should be considered when using this method to calculate a length on unprojected data?


texts: What functionality exists to calculate Cut and Fill for excavation?

texts: Does qgis have the ability to compare and quantify existing versus proposed elevations?


texts: leaflet count by featuretypes in featuregroup

texts:
Is there any way in leaflet to get the number of rectangles in a feature group?
I know this code that work: drawnItems.getLayers().length; but it counts all the objects within the drawnItems featuregroup. I only need the a specific type of features, for example rectangle


texts: Converting coordinates to something understandable for AutoCAD

texts:
I got this from a title and need to convert it and use it in AutoCAD with coordinates. Is there anyone who can tell me how to do this?

BEGINNING AT A POINT MARKED "1" ON PLAN,
BEING S. 29 DEG. 29’W., 1792.86 M. FROM BLLM NO. 1, ANTIPOLO CADASTRE; THENCE
THENCE S. 85 DEG. 56'E., 21.00 M. TO POINT 2;
THENCE S. 4-DEG. 04'W., 12.00 M. TO POINT 3;
THENCE N. 85 DEG. 56'W., 21.00 M. TO POINT 4;
THENCE N.4 DEG. 04'E., 12.00 M. TO THE POINT OF BEGINNING,


texts: Need to find pair-wise intersection for multiple polygon layers?

texts:
I have 20 different polygons and I'd like all the pairwise intersection among them.  I might be wrong, but I think Intersect isn't working because it's trying to find the intersection of all polygon layers?  But Pairwise Intersect will only accept two layers.  I've seen this similar question,

Is it possible to intersect multiple polygon shapefiles to identify areas where two or more layers overlap?

but manually counting the number of FID's in the union attribute table is not an efficient way.  I'm hoping there is a quicker option?

My attribute table looks like this:  -1, 0, -1, 0, 0, 1.  Where I'm assuming that would indicate the absence of the original two polygons and the presence of the third.


texts: How can I preserve date field format when importing to QGIS?

texts: I have a slight learning curve working here, is there a work around to importing spread sheet data without changing the type of data in the fields? i am importing some water test data along the rio grande and everything comes in fine except for the date field. date is formatted xx/xx/xxxx.


texts: GEE regression error: Image.reduce: ReduceBands reducer must have one input:

texts:
I'm trying to run a basic linear regression between two variables in Google Earth Engine. I have a two-band image: https://code.earthengine.google.com/?asset=users/stephenawood/reprexdata

To run a basic regression between the two variables and plot the results, I use the following code:

var linReg = reprexdata.reduce(ee.Reducer.linearFit());
Map.addLayer(linReg);

This generates the following error:

Layer 1: Layer error: Image.reduce: ReduceBands reducer must have one input: <Reducer>.

Note, I get the same error if I try print(linReg); rather than to map it. If I only execute the line starting var linReg I don't get any error.

I can use other reducers without error. The following works without an error:

var means = reprexdata.reduce(ee.Reducer.mean());
Map.addLayer(means);

This seems like it must be an issue of not being able to compute the regression reducer for some reason that I can't figure out.


texts: create spatial view in postgis on an ArcGIS table with ST_Geometry type

texts:
I am working on PostgreSQL 9.4.5 with Postgis 2.2.3.

I work on a spatial table created on ArcGIS, with geometry type "ST_POINT". I would like to create a spatial view of that table in postgreSQL using aggregation on few columns. However, I can't figure out how to write the query on how to aggregate my ST_POINT geometries into a multipoint geometry (ST_MultiPoint) in my spatial view.

Below is the script I am running, without success. The non-spatial aggregation is working perfectly fine for what I need (I tested it removing the geometry column) but I keep running into errors on whatever I try to write for aggregating the points into multipoint.

CREATE OR REPLACE VIEW global_ops AS
 SELECT string_agg(a.country, ', ') as country,
 a.rb,
 a.start_date,
 a.end_date,
 a.op_type,
 string_agg(a.activity, ', ') as activities,
 SUM(a.planned) as total_planned,
 a.links,
 a.year,
 a.project_num,
 SUM(a.actual_transfers) as total_actual,
 string_agg(a.transfer_modality, ', ') as modalities,
 sde.st_multipoint(st_AsText(''||a.shape||'')) as shape
 FROM global_ops_2018 a GROUP BY a.rb, a.start_date,
 a.end_date, a.op_type, a.links, a.year, a.project_num, a.shape ORDER BY
 a.project_num ASC;

This line here

sde.st_multipoint(st_AsText(''||a.shape||'')) as shape

is the one I'm trying to solve. Any ideas?


texts: Selecting connected lines from starting point using QGIS

texts:
I have a layer representing Watermain's made up of 2D lines. All the lines should be connected at the end points, but there are often very small gaps between them. Dead ends are connected at only one end but the vast majority are connected at either end. Sometimes a line may connect partway along another line, but this is incorrect as the ends should connect. Some lines may cross without connecting to reach other, and this is fine as long as the ends are connected.
Basically I want to be able to select all the lines connected to the starting line, and by inverting the selection identify those that need correction. Ultimately all pipes would be connected, and selected. I am not interested in the shortest route, but all the connections
I tried the topology checker in QGIS, but had trouble navigating to each pipe from the generated list. I would prefer a workflow of starting at a pipe and working progressively through the network until it was all connected.
Is there any commands to achieve this, or perhaps a plugin?


texts: Does anyone know of satellite imagery, such as Bing, that I can use to create my own polygons in QGIS by tracing?

texts:
For my study area, in western England I would like to effectively trace over certain fields to identify them as intense agricultural fields (such as bright green ones).

I know Google Earth and Google Maps forbid you to trace their maps so could I use Bing, OpenStreetMap, MapQuest layers that can be accessed by the OpenLayers plugin?


texts: Need to run a Groundwater Analysis Based on Points to a Raster

texts:
The title sounds daunting, but this task really isn't. I have well samples over a decade in certain locations through an area. These samples contain a value of that element at that point in time. I need to create a raster from these points to show the estimated dispersement in the area.

I know there is probably a lot of science behind this, but I need a quick and dirty way to do this. My background is much more in elevation in terms of raster, so I'm having an issue finding the right ArcGIS Tool for the job.

Here is the point's that I have, each colour representing a different range of values. It's also really important to note that only the latest sample is symbolized on the map when in fact there are four or five samples at each point that could have different values.

Like I said, I have a background in imagery / elevation when is comes to rasters. I needed a quick and dirty way to make this work, so I did a Natural Neighbor (Spatial Analyst) Point to raster analysis. I think it worked decently:

Really though I don't know. Does this tool fit the job or am I really stretching here?


texts: How to hide a layer dynamically from layers groups in Geoserver?

texts: I have created a layer group (using 15 layers) in Geoserver. I can access this layer group through Openlayers & it's working fine. As per the client requirement, I need to develop a layer switcher (for ON & OFF layer) for certain layers. So, How can I hide certain layers from this layer group. kindly help.