texts: Installing QGIS under Mint 18.2 (Ubuntu 16.04 Xenial)?

texts:
The problem

Following the instructions from the QGIS page, after the corresponding ppa channels are installed, the command

 sudo apt-get install qgis

fails with a message:

The following packages have unmet dependencies:
qgis : Depends: python-qgis (= 1:2.18.18+24xenial) but it is not going to be installed

Trying again with:

sudo apt-get install qgis python-qgis

leads to the error:

python-qgis : Depends: python-qscintilla2 but it is not going to be installed

Finally, trying with:

sudo apt-get install qgis python-qgis  python-qscintilla2

gives:

python-qscintilla2 : Depends: sip-api-11.3 but it is not installable

After some hours looking for a solution, I figured out that xenian only supports the sip-api-11.2. For some bizarre reason, QGIS packages point for the version 11.3, which does not exist in xenian, crashing the installation.

The Fix

Fortunately, a fix that worked on my machine was:

Download these individual .deb packages:

python3-pyqt5.qsci

python-qscintilla2_2.9.1+dfsg-4build1_amd64.deb

The package python3-pyqt5.qsci will force the installation of sip-api-11.2.

Install these packages in that same order:

sudo gdebi python3-pyqt5.qsci_2.9.1+dfsg-4build1_amd64.deb
sudo gdebi python-qscintilla2_2.9.1+dfsg-4build1_amd64.deb

After this qgis should install normally.

sudo apt-get install qgis


texts: Error 4 in rasterizing a polygon shapefile

texts:
I'm using QGIS Wien 2.8.1 trying to convert a polygon shapefile to a raster, but I get the following error in a qgis-bin pop-up:

ERROR 4: `' does not exist in the file system,
and is not recognised as a supported dataset name.

I have the 'gdal_rasterize' plug-in installed and I've tried executing this with the data and program files in the same directory (as well as multiple others) to make sure it's not a directory issue.

Any ideas?  I'll include the link to the polygon shapefile just in case it might actual be some issue in the data.
http://pubs.usgs.gov/of/2000/of00-444/of00-444_4c.tar.gz

Also, another question:  It seems that I can only rasterize upon numeric fields.  Is rasterizing upon string fields not possible?


texts: Corrupted exported image - Google Earth Engine

texts:
I am exporting a LANDSAT image for Malawi, but the exported image is corrupted, which I really do not understand why.

The link to the Code

https://code.earthengine.google.com/24db16ff4ecc2210f2f2c5e892204247

or here directly, this is my code for exporting the image

Export.image.toDrive({
  image: Malawi,
  description: '10_M',
  scale: 30,
  fileFormat: 'GeoTIFF',
  maxPixels: 120568920794 ,
  folder : '2010_Malawi'
});

and my code for constructing the image (just using the RGB bands)

var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw');
var Malawi = countries.filter(ee.Filter.eq('Country', 'Malawi'));
var landsat = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');

var Malawi = landsat
.filterBounds(Malawi)
.filterDate('2010-01-01', '2010-12-31')
.select(['B1', 'B2', 'B3'])
.median();


texts: Why do file geodatabase and SDE readers produce different fme_geometry types?

texts:
I have one very simple workbench that I am running against both SDE (9.3) and FGDB (10.0) sources.

I noticed two different fme_geometry values for same type of objects.

For the first one, i.e SDE 9.3, I have fme_geometry=fme_donut, but for the same features in FGDB 10.0 I have fme_geometry=fme_aggregate.

Is this a "problem" of converting ESRI data from 9.3 to 10.0 or how FME interprets geometry?

Because in case of fme_aggregate, FME is doing "additional" job of constructing donuts again:

Calculating inside points
Determining uppermost containment
Constructing donuts from xxx points

(Visible in log window)


texts: Getting EPS vector as georeferenced layer into QGIS

texts:
I have an eps-file with polygons of landuse. The file has coordinates to georeference it. I need to convert it into a geo vector.

I could do the following steps:

export/convert the eps to raster (unsing Inkscape)
georeference it with GDAL or J Rectifier (great tool!)
reclassify it (using r.reclass from GRASS) and
vectorize it with GDAL, GRASS or SAGA.

but is there another way without rasterizing the eps-vector?


texts: Fnding wildcard string using Field Calculator script code?

texts:
I'm trying to find part of a text string in the CLASS_NAME field and return a text string to the Abrv_Name field using the Field Calculator in ArcGIS 10 via a Python parser.

Here is what I have tried among other things with no success.

This worked within the field calculator.


texts:
I would like to speed up the process corresponding to this command:

gdalwarp -t_srs EPSG:4326 -overwrite input.ntf output.tif

Is there a way to use parallel processing in GDAL?
Semi-offtopic: If not, do you recommend a non-GDAL solution to speeding up a gdalwarp-like process?

Here are the web pages I've looked at:

Is it possible to do parallel processing in GDAL and QGIS?
[gdal-dev] GDAL raster processing: parallel computing

Edit: I asked this question because I thought I saw my CPU% stay below 100% when processing gdalwarp. However, on a second look, it reached 555.5%.


texts: Implementing a classification on Google Earth Engine

texts:
I am attempting to implement a classification based from a decision tree (attached as image below) where there are 5 classes. These classes are based off the different water reflectance values in Landsat 8 bands (namely NIR, Red and Green). I have no idea how I should go about doing this. Someone suggested using image.expression however it has not worked for me. This is the code I've used thus far:
(landsat is landsat 8 surface reflection)

var image = ee.Image(landsat.filterDate ('2018-01-01', '2018-03-01').filterBounds(point).sort('CLOUD_COVER').first());

//print image
print('a landsat scene', image);

//implement MNDWI

var green = image.select('B3');
var swi = image.select('B6');
var mndvi = green.subtract(swi).divide(green.add(swi)).rename('MNDWI');

var mndvipara = {min: 0, max: 0.7, palette: ['white', 'blue']};

Map.addLayer(mndvi, mndvipara, 'MNDVI image');

//classification attempt
var lpb =image.expression (
'nir > 0.018', {
'nir': image.select ('B5')});
Map.addLayer (lpb, {min: -1, max: 1, palette: ['FF0000', '00FF00']});

[1]: https://i.stack.imgur.com/TKSbZ.png (if image doesn't work, classification can be found at https://www.sciencedirect.com/science/article/pii/S0303243415000641 - Figure A.1


texts: Find Minimum Value of Point within a buffer output as a point file

texts:
I have a series of 30m buffers along a polyline (river). Within each buffer I have a series of points each point represents a 10m x 10m cell with an elevation value known as "Elev".

I want to find the one point with the minimum elev value within each buffer and output it as a point file. In a future geoprocessing set I'll like to snap this point back to the river.

ArcGIS 10.1, arcinfo license.


texts:
I have 273 features in the feature class and i want to select them one by one and use the statement in the following code.

At the line:

*arcpy.SelectLayerByAttribute_management (fc1,"NEW_SELECTION","OBJECTID=79")*,

The feature whose object id is 79 one of the features. When i use this code, I do it for only that feature. I want to apply the statement to all features one by one with a loop. I try to do it, but i am new for programming with python.

import arcpy
from arcpy import env
env.workspace = "D:\pomme de terre 2\pomme de terre2.gdb"
env.overwriteOutput = 1

fc1 = "Flowlin24K_sl_Uns"
fc2 = "New_Poyline"

arcpy.SelectLayerByAttribute_management (fc1,"NEW_SELECTION","OBJECTID=79")
arcpy.SelectLayerByLocation_management(fc1,"INTERSECT",fc1)
arcpy.SelectLayerByAttribute_management (fc1,"REMOVE_FROM_SELECTION","StreamLevel>6")
arcpy.CopyFeatures_management(fc1,"hidden")
arcpy.SelectLayerByAttribute_management (fc1,"CLEAR_SELECTION")

coordList = []

arcpy.FeatureVerticesToPoints_management("hidden", "vertices", "MID")
arcpy.AddXY_management("vertices")
rows = arcpy.SearchCursor("vertices")

for row in rows:
    X = row.getValue("POINT_X")
    Y = row.getValue("POINT_Y")
    coordList.append([X, Y])

del row, rows

coordList.sort()

point = arcpy.Point()
array = arcpy.Array()

for feature in coordList:
    point.X = feature[0]
    point.Y = feature[1]
    array.add(point)

polyline = arcpy.Polyline(array)

array.removeAll()
arcpy.Delete_management("vertices")
arcpy.Delete_management("hidden")
arcpy.CopyFeatures_management(polyline, fc2)


texts:
I'm using QGIS 3.4.  I have a point shapefile representing population density, and I would like to identify the 12 most optimal locations for service centers that would  serve the maximum amount of people in a given area. (See screenshot)

I tried the "K means clustering tool", which essentially plugs through multiple iterations that assign each point to the nearest cluster until the optimum number of clusters requested (in this case 12) are reached; However, the tool doesn't seem to be working, the outpoint is a duplicate of the input points shapefile.

I've used Network Analyst for a similar project in ArcGIS (see screenshot, the firework clusters are the optimal locations identified as a result of the cluster analysis), but I don't know how to proceed in QGIS.  Suggestions?


texts: Create a subset of geometries which intersects with each other

texts:
I have a set of geometry in my PostgreSQL table that consists of Points, Linestrings & Polygons. From this set I need to create a subset of geometries which intersects with each other.
I need a method without looping through each one of them as there can be more than 1000 geometries at once.
I need to fetch the geometries that are marked in red in this image


texts: Alternative of ST_Translate function in SQL SERVER 2016

texts:
I am converting a PostgreSQL script into SQL SERVER 2016.

In PostGIS we use ST_Translate as given below :

geometry ST_Translate(geometry g1, float deltax, float deltay);

Is there any alternative of ST_Translate PostGIS function in SQL SERVER 2016 ?


texts:
I have a data set at administrative level (generated by appending csv to shapefile of admin areas) that i want to clip to physical building features (from Openstreetmap). To tackle this i've approached as below

Created a layer file, turned this into a representation then used the Analysis clip tool on that file with the building out lines and reapply the symbology from the original file to retain the correct colours. The problem i'm getting is that the output appears incorrect below is the overlay of original data in colour with building outlines in black;

Using the method described this is what Arcmap plots, it appears to be doing something with the data as it plottting different colours on the clipped outlines to the underlying values shown on the original unclipped areas.


texts:
In qgis 3, postgis, no DEM and zero python knowledge.
I've got linestring 3d but Z coordinates is empty, Z-values are stored in table.
Z1 for the first node, Z2 for the last node.
I found a lot of ways to extract Z to populate the table but never opposite.
I try "set Z value" processing, but I can't assign value Z1 to first point and Z2 to last point.
I turn around with start_point, end_point,st_z, maybe st_makeline and others but I can't do it.


texts: Using Quickbird to discriminate vegetation response and spectral signatures - which season?

texts:
I have a study area in N. Florida loblolly pine plantation with 3 treatments (drought, fertilization, drought+fertilization, control) and 4 replicate plots. All plots will be in the same scene. These treatments are expected to alter the growth and response of pine trees.

The only available existing imagery at the initial time of treatment installation is from quickbird (one in March and one in December) - this will be the baseline image. Ideally, it would be nice to acquire imagery every year and season, it's not feasible due to budget constraints. So, I will probably purchase one, maybe two more images over the next 1-2 years. They should be the same season as the "baseline" scene.

My Question:
Which season, if any, is best to discriminate spectral signatures between treatments/plots to really tease out differences in growth due to treatments and estimate NDVI, LAI, and other vegetation indices? I would suspect winter, when everything is not as lush and green and somewhat stressed (relative to summer). What is your expert opinion? Thanks in advance.


texts: Unable to import shapefile data (both POLYGON and MULTIPOLYGON) to postgis using GDAL osgeo, ogr

texts:
I am trying to import shapefile data (polygons with all of the user fields) to postgis using osgeo. My code fails when shapefile contains both POLYGON and MULTIPOLYGON data.

from osgeo import ogr, osr

# Open shapefile
ogr_sf = ogr.Open('shapefiles/tst_2.shp')
shapefile_layer = ogr_sf.GetLayer(0)

# Open database connection
ogr_db = ogr.Open("PG:dbname='geo_development'"
                  "host='postgres' port='5432'"
                  "user='user'"
                  "password='password'")

db_layer = ogr_db.CopyLayer(shapefile_layer, 'tst_table['OVERWRITE=YES'])

I've got the next error message:

  Warning 1: Geometry to be inserted is of type Multi Polygon, whereas the layer geometry type is Polygon.
  Insertion is likely to fail
  ERROR 1: COPY statement failed.
  ERROR:  Geometry type (MultiPolygon) does not match column type (Polygon)
  CONTEXT:  COPY tst_table, line 2, column wkb_geometry: "0106000020E6100000020000000103000000010000000700000083F9A3C6597B5EC0ACADB807A7AE424040B2E29AEC7A5EC0..."

I've tried to walk around defining the postgis table first and then copy data:

from osgeo import ogr, osr

ogr_sf = ogr.Open('shapefiles/tst_2.shp')
shapefile_layer = ogr_sf.GetLayer(0)
layer_definition = shapefile_layer.GetLayerDefn()

ogr_db = ogr.Open("PG:dbname='geo_development'"
                  "host='postgres' port='5432'"
                  "user='user'"
                  "password='password'")

# Define spatial reference system
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# Create layer
db_layer = ogr_db.CreateLayer('tst_table', srs,
                              ogr.wkbMultiPolygon, ['OVERWRITE=YES'])
for i in range(layer_definition.GetFieldCount()):
    db_layer.CreateField(layer_definition.GetFieldDefn(i))

# Copy features
for feature in shapefile_layer:
    db_layer.CreateFeature(feature)

This time the error msg is opposite to the previous one:

  Warning 1: Geometry to be inserted is of type Polygon, whereas the layer geometry type is Multi Polygon.
  Insertion is likely to fail
  ERROR 1: COPY statement failed.
  ERROR:  Geometry type (Polygon) does not match column type (MultiPolygon)
  CONTEXT:  COPY tst_table, line 1, column wkb_geometry: "0103000020E6100000010000000C000000116AC8BB107B5EC043A499E388B04240D980B833AC7A5EC07653090DC6B04240D8..."

Does anyone know how to overcome this?


texts:
What are the rules/guidelines for what the total population in a US census block group has to be to justify a split?

I'm wondering if the US Census Department has any guidelines about when to split a census block group into separate parts.


texts: How to choose a GRASS kernel density function?

texts:
I'm starting to use the Grass module on QGIS, and I want to do a heatmap, so, I'm creating a new surface using the kernel function. But, GRASS come with 7 types of kernel functions, how do I know which one to choose?

Regards


texts: How best to relate historical testing data to a shapefile?

texts:
I am tasked with bringing what has historically been only tracked in Excel into a shapefile. I am using ArcMap, but also have QGIS installed on my computer

What these shapefiles will include is inventory and location information about the road segment being tested which is easy enough to do.
The problem I have is that the agency I am doing this workfor also wants all their historical testing data to be included in these shapefiles (done yearly and there are about 10 numbers to keep track of each year). Now the only thing that I can think of is adding a prefix to indicate the year (2002_Friction or something) and then each column of testing, but that would result in a ridiculously huge number of columns. I am going to have to probably do this several dozen times at the minimum so I'm trying to find an efficient way to do this.


texts: Joining X and Y coordinates form CSV with shapefile

texts:
I am working in QGIS. I have a shapefile showing regions of Sweden.
In a .csv-file I have X and Y coordinates which I want to link to the shapefile spatially. How can I do it?


texts: Creating Swath View from Point data using Python?

texts:
My situation

Input: point data from planting file

Feature Attribute:
    [('DeletionFlag', 'C', 1, 0), ['Field', 'C', 254, 0], ['Dataset', 'C', 254, 0], ['Product', 'C', 254, 0], ['Obj__Id', 'N', 15, 4], ['Time', 'D', 8, 0], ['Duration_s', 'N', 15, 4], ['Area_Count', 'C', 254, 0], ['Elevation_', 'N', 15, 4], ['Swth_Wdth_', 'N', 15, 4], ['Y_Offset_f', 'N', 15, 4], ['Distance_f', 'N', 15, 4], ['Track_deg_', 'N', 15, 4], ['Rt_Apd_Ct_', 'N', 15, 4], ['Tgt_Rate_k', 'N', 15, 4], ['Pass_Num', 'N', 15, 4], ['Seed_Cnt__', 'N', 15, 4], ['Plant_Pop_', 'N', 15, 4], ['Speed_mph_', 'N', 15, 4], ['Prod_ac_h_', 'N', 15, 4], ['Prdt_Amt', 'N', 15, 4], ['Date', 'D', 8, 0], ['Population', 'N', 15, 4], ['Product___', 'C', 254, 0]]

Sample data

https://drive.google.com/drive/folders/1RLflZNFrN0n0Ew35hTTyorecvo-l_ISq?usp=sharing

I have successfully created polygon from point data using below link
Code:
Creating Polygon Shapefile from Point data using QGIS?

import shapefile
import math
fileName='planting';
sf = shapefile.Reader("planting/" + fileName )
print(sf.fields)
w = shapefile.Writer('swath/' + fileName)
w.field('Product', 'C')
def rotate(point, angle):
    angle *= math.pi / 180.
    px, py = point
    qx = math.cos(angle) * px - math.sin(angle) * py
    qy = math.sin(angle) * px + math.cos(angle) * py
    return qx, qy

def translate(point, translation):
    return [point[0] + translation[0], point[1] + translation[1]]

def rectangle(width, height, angle, center):
    p1 = (-width/2, -height/2)
    p2 = (-width/2, +height/2)
    p3 = (+width/2, +height/2)
    p4 = (+width/2, -height/2)
    rotated = [rotate(p, angle) for p in [p1, p2, p3, p4]]
    translated = [translate(p, center) for p in rotated]
    return translated

for feature in  sf.shapeRecords():
    width = 0.8 * 0.3048 * feature.record['Swth_Wdth_'] / 100000  # MAGIC
    length = 1.2 * 0.3048 * feature.record['Distance_f'] / 100000  # MAGIC
    direction = feature.record['Track_deg_']
    rate = feature.record['Rt_Apd_Ct_']
    product= feature.record['Product']
    center = feature.shape.points[0];
    points = rectangle(width, length, -direction, center)
    w.poly([points])
    w.record(product)

Point data

Swath view

this is implemented in python

I need to find

The swath view we have achieved from the point data in the planting file is not as desired. On certain areas in the field, some sensors in the swath are turned off, hence the effective swath width decreases. In our case, this is being plotted on top of the existing swath row.

How can we handle the varying swath widths when sensors are turned off?


texts: Calculating coordinates for halfway around world from given location

texts:
I am curious about how to mathematically determine "half-way around the world" from a given spot.
For example, from where I live in Monterey, what is half-way around the world heading either due east or due west?
I also want to know what is halfway around the earth the other direction - heading due north or south.
The coordinates for Monterey are approximately latitude 36.6 and longitude -121.89
My idea is that halfway around the earth from Monterey, heading either east or west, is 36.6 latitude (staying on the same plane/parallel of latitude) and 58.11 longitude, in other words in Birjand, Iran. I'm pretty sure I'm right about that.
But I'm having more trouble with the other direction. Is latitude -53.4 and longitude -121.89 (same meridian of longitude as Monterey) correct, and if so, where is that, exactly?
When I enter those coordinates in both Bing Maps and Google Maps, it just shows me a spot in the ocean, but won't allow me to zoom/expand out from there so as to give me my bearings (no pun intended) or context.
What would be the mathematical formula to apply here?


texts: To write Error message grouped by a common field attribute value and check another field value is matching for that same group instance using arcpy

texts:
The name of my layer is "spat_join_output". So I have a field named 'label' and another two fields 'end_point_x' and 'end_point_y' for the same label values i want to check if the end_point_x and end_point_y values are same , i want to drop an error message in another field named as "Error".

I have my code as :-
array_fibre = arcpy.da.FeatureClassToNumPyArray("spat_join_output", ["label", "end_point_x", "end_point_y", "OBJECTID"], skip_nulls = True)
        # arcpy.AddMessage(array_fibre)
        for data in array_fibre:
            arcpy.AddMessage(str(data))

    if arcpy.Exists("spat_join_output"):
        arcpy.AddField_management("spat_join_output", "Error", "Text", 100)

    with arcpy.da.UpdateCursor("spat_join_output", ["label", "end_point_x", "end_point_y", "Error"]) as cursor:
        for row in cursor:
            if row[0] != None and row[1] != None and row[2] != None:
                data = self.find_endpoints(array_fibre, row[0])
                if len(data) > 0:
                    arcpy.AddMessage(data[0])

        del cursor

I am trying to build a function for this, using an array method,
def find_endpoints(self, array, id):
    count = 0
    array_id = []
    end_point_X_y = []
    array_label = []
    if len(array) > 0:
        array_count = 0
        for data in array:
            if len(array_label) == 0 and id in data[0]:
                array_label.append(data[0])
            elif len(array_label) > 0 and id in data[0]:
                array_label.append(data[0])

This is how my array looks:-
('CAS/C003', '-877644.7019000016', '6893765.545000002', 2)
('CAS/C003', '-877939.7175000012', '6893358.312799998', 3)
('CAS/C003', '-877244.6374999993', '6893426.772500001', 4)


texts:
I am working on calculating a slope raster in Mapinfo using a 15cm DEM derived from a LiDAR point cloud and the output data will have 3 predefined breaks displayed so I need to blur those relatively small areas across a larger DEM area fairly accurately.

I've read about tools which use FFT or Gaussian Blur but can't find anything for Mapinfo (64 bit). Can anyone suggest a tool or command that can be used?

edit from comments:

I need to modify the 50cm pixels from their existing value to a value based on a number of surrounding pixels to generalise the slope. I believe this could be achieved using Fast Furier Transform or a Gaussian blur algorithm but I am not aware of which tools/commands to use in Mapinfo 64bit to accomplish this.


texts: Determining spatial cluster of park types in ArcGIS Desktop

texts:
Question 1:
I am trying to determine which (if any) of my different park types are spatially clustered across my study area (city boundary).

There are 162 total parks, 9 nominal categories (which represent statistically determined clusters or park 'types' based on 33 social, built, and ecological variables).

Also, would I run this analysis on the polygons (park boundaries) or points (based on centroids), or some other way?

Question 2 (related but different):
I am also wondering if I can run the entire analysis (including spatial) in Arc? The cluster groups I have now reflect the results of an SPSS process including
1. a factor analysis on the 33 variables to reduce redundancy then
2. a two-step cluster analysis (Hierarchical to determine clusters then k-means to form groups).


texts: MapServer: how to publish a TMS using WMS?

texts: Is it possible with MapServer to publish a WMS from a local TMS generated by gdal2tile?


texts: google earth profile used in mashup or flex

texts:
I really like the way google earth profile works,
 1. mouse over line, right click, select profile
 2. slide the mouse over the created 3d profile an dsee where that is on the map.

Has anyone an idea of whether this capability could be mashed outside of GE (I probably already know the answer)? Or if there is anythig esle as slick that could be integrated with flex or arcgis server?

Thanks


texts: Select maximum values from list of raster layer : ArcMap10

texts:
I have 10 raster layers. Each layer has the dimension of 100 columns x 100 rows and actually are temperature values from year 2000 to 2010.

I want to ask ArcMap to produce a raster layer that represent the maximum temperature value from 2000 to 2010.

Is there any function in ArcMap10 that suite to do that?


texts:
If while I'm executing gdalinfo against a raster layer I get a positive value for pixel size Y, what of these scenarios is possible?

The raster's map coordinate system sets the origin of the data in the upper left corner.
The raster is not georeferenced and its coordinates are interpreted as image pixels, instead of map units.
Any other situation (corrupt file?, other?)

I know this value is almost always negative. For example, in case of UTM coordinate system. But, how common is to have a positive value for Y coordinate?

Many thanks in advance


texts: Creating ArcGIS Service (FeatureAccess) with dynamic datasource?

texts:
I am using ArcGIS Desktop and ArcGIS Server version 10.3. I have also created an enterprise geodatabase in Oracle 12C.

In the SDE schema I have a table containing a road network (about 10000 links), for the sake of brevity let's assume that the table has the following structure:

LINKS_TABLE

ID (NUMBER(10,0)), SHAPE (ST_GEOMETRY)

Furthermore, I have a second table which holds link-related attributes (e.g. vehicles per hour). In this second table I store data for different scenarios (e.g. vehicles per hour in weekdays, in weekends, during summer, during winter, etc). The table structure is the following:

SCENARIO_TABLE

ID (NUMBER(10,0)), LINK_ID (NUMBER(10,0)), SCENARIO_ID(NUMBER(4,0)), SCENARIO_VALUE(NUMBER(10,0)).

I want to display the aforementioned network on a web application using either leaflet (ESRI Leaflet) or OpenLayers. I need to be able to select different scenarios (most likely more than one at the same time) and display them as different layers. This means the datasource for each layer will be something like

SELECT LINKS_TABLE.ID, LINKS_TABLE.SHAPE, SCENARIO_TABLE.SCENARIO_VALUE FROM LINKS_TABLE INNER JOIN SCENARIO_TABLE ON LINKS_TABLE.ID = SCENARIO_TABLE.LINK_ID WHERE SCENARIO_TABLE.LINK_ID = XXX

I want to deliver the output as an ArcGIS FeatureAccess layer.

The SCENARIO_TABLE will constantly be updated by users (or by background processes), so I cannot export the required services beforehand.

How can I create an ArcGIS Service where I can pass the scenario_id as an argument during runtime?

The web application will most probably (99%) be written in ASP.NET using C#.


texts: Building Socket Server to receive GPS location data

texts:
I am working on a web application which receives data sent from trackers via TCP Packets containing location information. Up until this point we have been working with dummy data within the mysql database but we are now ready to start building the server which will receive data from the remote devices.

As we are just exploring our options at the moment, I thought I would see does anyone have any suggestions on the best route to take when setting up a server to receive TCP packets of data and then store them in a remote database?

Our general requirements are as follows:

Hosted on Amazon EC2
Allow multiple concurrent connections to receive data (potentially in the thousands)
Store data in remote MySQL database


texts: Counting polygons in one layer that are overlapping/touching/contained in each polygon in separate layer of contiguous polygons using QGIS?

texts:
I have a layer of contiguous (pink) polygons (UK Parliamentary constitutencies) and a separate layer containing (blue) polygons representing protected nature areas (Natura 2000 areas):

Each constituency consists of a single polygon. Each Natura 2000 area can consist of one or more polygons.

For each constituency I need a combined total for the number of Natura 2000 areas that are:

contained in
touched by
overlapping

that constituency.

I've added both as vector layers in QGIS and tried using Spatial Query and Topology Checker plugins and Geoprocessing tools (Clip, Intersect). But I don't really know what I'm doing. Although I was able to select constituencies where there is an intersect I couldn't see the number of intersects per constituency and this didn't solve my need to include touching and contained polygons in the totals.

And just to be clear about Natura 2000 areas consisting of one or more polygons. In this screenshot (contained in a box) is an example of a single Natura 2000 site consisting of multiple polygons. Because all of these polygons are contained in a single constituency we'd only want to count them once as a single Natura 2000 area. But if they were spread across the boundaries of constituencies we'd count them for each constituency any part of them touched/overlapped/was contained by. In the same way we'd count the single polygon Natura 2000 site that overlaps two constituencies (circled) twice - once for each constituency it covers.

(For background, I'm campaigning to protect Natura 2000 areas and my aim is to develop a target list of UK Members of Parliament prioritised by those who have the highest number of these areas in their constituency.


texts: How to use the Cartographic Composer of GRASS 6.4.3?

texts:
Following on from my previous question (which was, again, how to use the Cartographic Composer)...

I've got a GIS in GRASS. I need to print various bits of it. Problem is, it needs images on it whenever we print it - the logos of our clients and stuff.

Now, I've managed to get the Cartographic Composer to show the logo images, now I just need it to show the map - or more rather, the Map Display and what that is currently showing. As everyone who will be using the GIS is used to ArcMap (but it would cost too much to get the licences), we basically want it so that the user can just press the print button and GRASS will print out the Map View, plus the logos.


texts:
The ultimate question, I need a transformation matrix to take a point in local space representing a roughly 500m x 500m place in New Mexico centered at -108.619987456 long 36.234600064 lat.  The final output needs to be in geocentric coordinates and I can during initialization get any number of points on the map such as the corners, center, western most along the center etc.  During run time I would like to have a transformation matrix that can be applied to a position in local space expressed in meters and get back an approximate GCC coordinate.

The center of the point in GCC: -1645380 -4885138 3752889
The bottom left corner of the map in GCC: -1644552, -4881054, 3749220

During run time I need to be able to multiply (-250, -250, 0) with the transformation matrix and get something close to (-1644552, -4881054, 3749220).

I've spent more than a week trying to research a solution for this and so far nothing.

Given an identity matrix as our starting position and orientation. Then using geotrans and the known lat, long, and height of the starting position I get an x,y,z geocentric coordinate. A vector from the origin of (0,0,0) gives both the up and translation for the matrix. However, I need the forward and right so that I can pass a distance in meters from the origin into the transformation matrix and get a roughly accurate GCC. Do I have all of the inputs I need to calculate the right and forward?

Inputs
Origin: 0, 0, 0
Global: -1645380, -4885138, 3752889
Up (Normalized): Global - Origin

Desired Outputs
Right: ? ? ?
Forward: ? ? ?

So with the right and forward added to the up and translation I already calculated I would have a transformation matrix.  I could then apply the matrix to a vector of say (50, 50, 0) and get something within about 0-3 cm of the output if I feed the lat/long back into geotrans.  This matrix would only ever be used small maps of about 500m x 500m so the curvature of the earth is negligible.

In reply to whuber,

I don't know your level of experience so I will start with the very basics.

An identity matrix is a 4x4 matrix that you can multiply by a vector and it returns the vector unchanged. Such as below.

1 0 0 x=0
0 1 0 y=0
0 0 1 z=0
0 0 0 1

The x, y, and z of the matrix are how much to translate the vector by and the first 3 rows and columns are the rotation.  Below is a tranformation matrix that does nothing to the orientation of a vector but does translate.  This is what you would want for two axis aligned worlds in different locations.

1 0 0 13
0 1 0 50
0 0 1 -7
0 0 0 1

If you were to multiply a vector of (10, 10, 10) with the above transformation matrix you would get an output of (23, 60, 3).

My problem is the axes are not aligned and the "plane" of the world I am trying to get the local coordinate converted to is projected on the ellipsoid of the earth.

Geotrans is library that you can use to pass a coordinate from one system such as geodetic or GDC (lat, long, height) and get back another such as geocentric or GCC (x,y,z).

For example:
If my game world was representing a map centered exactly on the prime meridian and equator the transformation matrix would look something like below. I am just guesstimating the Y and Z rotations and might have them flipped and/or the wrong sign.

0 0 1 6378137
0 1 0 0
1 0 0 0
0 0 0 1


texts: Converting polylines as paths to polygons in one shapefile using ArcMap

texts:
I have acquired raw path GPS data where seabird colony boundaries were physically recorded and uploaded in the form of polylines.
A lot of these polylines are a bit messy and overlapping, but there are about 10 non-contiguous colonies on the map area.
I wish to convert these into polygons under one shapefile to represent seabird colonies. I have attempted several different methods, but in converting polylines to polygons under one feature class file, I have a hard time with most of the polylines converting as if they are not perfectly connected the command will not recognize the polygon I wish to create.
Is there any easier way of smoothing and connecting polylines other than manually using the create features tool?


texts: Geoserver SQL View that uses an existing fuction

texts:
I am attempting to create a Geoserver SQL View layer that uses an existing function in my Postgres/Postgis database.

Starting with the simplest case I have successfully created a layer using the following SQL statement on the Edit SQL View page:

SELECT geom FROM map.parcel WHERE pkid = %pkid%

Then I create a function in the database:

CREATE OR REPLACE FUNCTION map.gs_parcel (_pkid integer)
RETURNS TABLE (geom public.geometry) AS
$body$
BEGIN

RETURN QUERY
SELECT p.geom FROM map.parcel p WHERE pkid = _pkid;

END;
$body$
LANGUAGE 'plpgsql'
VOLATILE
CALLED ON NULL INPUT
SECURITY INVOKER
COST 100 ROWS 1000;

And I try to create a SQL View using:

SELECT geom FROM map.gs_parcel(%pkid%)

But I get the error:

ERROR: function map.gs_parcel(integer) does not exist Hint: No function matches the given name and argument types. You might need to add explicit type casts. Position: 34

What am I doing wrong? I need to expand this simple example to take multiple parameters and add logic, which I can only do using pl/pgsql.


texts: Calculate sine and cosine in Google Earth Engine

texts:
I need to use the sine and cosine of an angle. How can I calculate it?

I need something like this:

var a = 40 (angle)
var sina = sin(a) -> in order to have in output a number


texts: Adding tracks to gps with QGIS feature attributes as name

texts:
I am attempting to send a vector layer to a garmin gps in QGIS. I am using the GPS Tools extension and have created a new GPX layer. I can add waypoints without any issues, but when I copy features from a line layer onto the track gpx layer, it will not save.

I can save this layer in a gpx format and open it with MapSource, but I'd like to set one of the attribute as the track name automatically as I have over a thousand.

Any help would be helpful. Also, if anyone could explain why the gpx layer won't save when I paste features into it that would be great.


texts: Getting IndexError: list index out of range from arcpy.mapping.ListLayers() in Python Window of ArcMap?

texts:
i want to turn on label using arcpy, i use current mxd and my mxd path in my script and my layers also active in the current MXD, but it give me this eror, can somebody help me?

here the script i used :

import arcpy

... mxd = arcpy.mapping.MapDocument(r"C:\Users\P51\Desktop\testingLABEL.mxd")

... df = arcpy.mapping.ListDataFrames(mxd,"Layers")[0]

... layer = arcpy.mapping.ListLayers(mxd, "point10", df)[0]

... layer.showLabels = True

and this is the eror :

Runtime error
Traceback (most recent call last):
  File "", line 4, in
IndexError: list index out of range

here is the screenshot :


texts: SAGA Tool in QGIS Error: Load table: Proj.4-WKT Dictionary... failed

texts:
I have some problems with running SAGA tools in QGIS. I've downloaded and installed the software via the OSGeo4W Installer. It downloaded QGIS 2.12 and Saga 2.1.2. I wanted to use the Tool "Shapes To Grid" to convert a Point Feature to a Raster but I've always got the Error that the "gdal_MrSID.dll" wasn't able to load although the dll was in the folder. So I tried installing the latest SAGA Version 2.2.2 64bit, set the right path in QGIS and then it worked!
Afterwards I tried to use the "Change Grid Values" Tool but it didn't work. I got some errors:

Load table: Proj.4-WKT Dictionary...
failed

And later:

Error: select a library

That is very strange because the Tool works perfectly with the path set to Saga 2.1.2. I guess it has something to do with path settings but I don't know how I could solve this.

Has anyone an idea what could cause this? I am very new to SAGA and QGIS. I used ArcGIS before. If you need any more information about the problem, my system or my paths I will tell them to you. Right know I don't really know what information could be useful for this topic.


texts:
If you want an accuracy of seven digits, should you then round to that value or simply just cut the numbers after that?
Example:

51.1936137567  to 51.1936138
or
51.1936137567  to 51.1936137


texts: Creating a continuous raster road network using Polygon to Raster in ArcGIS Pro?

texts:
I have a polygon forest road network that I would like to turn into a raster. I have attempted to use the ArcGIS Polygon to Raster tool, but whenever the majority of a cell falls outside of the polygon, it is given a value of No Data, meaning that the resulting raster does not have continuous roads. See image below, with black outlines representing the vector features.

I have tried setting a priority field, changing the cell assignment method, and buffering the roads (this works but produces inaccurate results). The only solution I have found so far is to calculate the center lines of the roads and then use the Polyline to Raster tool to create the road network. This works (see below), though it adds complexity to the workflow.

Does anyone know of a method of creating continuous road networks from polygons, without needing to calculate the center lines?

I am using ArcGIS Pro 2.2.3.

I have attempted to use a finer cell size (1m) and then re-sample back to the original cell size (10m). Unfortunately, while the output is better, the original problem persists (see below).


texts: Create ESRI ArcSDE Connection File (.sde) using Python/ArcPy at the ArcView license level?

texts:
Is there an alternative way to create ArcSDE connections in ArcGIS 10.0 using Python at the ArcView license level? I need to create 4 separate connection files for a MXD SDE data source re-point tool.

ArcView licenses are the only available.

The SDE connection files need to be created on-the-fly for each specific user so they can enter their own username and password.

The arcpy.CreateArcSDEConnectionFile_management method is only supported at the ArcEditor/ArcInfo license level.


texts: Arcpy - Loop through a feature class and create Euclidean Distance for every feature

texts:
I am trying to use SearchCursor to loop through a feature class and create a Euclidean Distance raster dataset for every feature in that feature class. This ED raster dataset will be used for further analysis. However, I got this error message which I cannot solve. Can anyone explain or help please.

cursor = arcpy.SearchCursor(Tank_Location_Data)

for row in cursor:

    arcpy.gp.EucDistance_sa(row, "50", "1", "#")

Error message:

arcpy.gp.EucDistance_sa(row, "50", "1", "#")
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\geoprocessing\_base.py", line 498, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool


texts: Preserve table info when converting AutoCAD .DWG to MapInfo .TAB

texts:
I've been using Universal Translator to convert some .DWG files into .TAB.
I would like to know if there is there is a way to preserve the data in the AutoCAD file's table when converting to .TAB?
Currently the tables are empty after the conversion:


texts: ArcGIS projected area greater than geodesic area

texts:
Question

Shouldn't the geodesic area always be greater than the planar area of a polygon?

Background

Tested in ArcMap and ArcGIS Pro, but several my polygons in a shapefile are having projected areas that are greater than the geodesic area.

I've loaded in a shapefile that has a GCS of WGS1984; the dataframe is set to GCS: WGS1984. Once loaded I perform the following:

Geodesic Area

Add Field to attribute table and set to double
Right click field and selected Field Calculator
Supply the following Python input: !Shape!.getArea("GEODESIC", "SQUAREMETERS")

Projected Area

Select the Project tool in Data Management > Projections and Transformations
Project the layer to PCS: World Sinusoidal

Then I use Python to check the 2 fields against one another. A quick check of:

# pseudocode
for each polygon:
    diff = geodesic_value - projected_value

A negative value for diff means that projected was larger; postive is the inverse. It was my understanding that geodesic would (and should) always be larger.


texts:
I'm a bit confused on how the zonal statistics tools work in QGIS (either the 'Zonal Statistics' tool under the Raster tools from the QGIS geoalgorithms toolbox or the 'Raster statistics from polygons' tool under the SAGA Vector<->Raster tools).

Basically I have a raster map layer and some small polygons in a separate layer. I want to assign each polygon the value of the raster below. In case a polygon overlaps multiple raster cells, I just want to have the maximum.

Yet the output doesn't make any sense and I can't get behind the reasoning of those results (which are obviously not the maximum of the underlain raster cell):

Any ideas what I'm doing wrong or what I didn't consider yet? The results with the python API in qgis.analysis.QgsZonalStatistics() (according to this question here) are much better, but it does only offer count, mean and sum...but for some reason no max or min.


texts: How to display geoJSON from an URL in a layer?

texts: I am looking for some type of plugin or code which would allow me to add a geojson URL, define the areas to input the lat/long variables from the visible Qgis canvas map and call/render the results within Qgis. Anyone know how to do this?


texts:
I've got a dataset of buildings roof segments in shapefiles. I need to sequentially number the roof segments from 1 to n of each building. Is there a way to do it with the QGIS field calculator?
Actual:

Expected:

It's pretty similar to that question on arcpy: Auto increment by group in ArcGIS ModelBuilder


texts: How to implement tile caching a WMS vector layer and use that cached tiles with leaflet?

texts:
I have created a vector layer using GeoServer. Now, i have a GeoServer WMS for that layer. I would like to implement tile caching. I tried TileCache tool. But i have encountered a problem then i asked the solution in this question, but nobody answered it.

When i use the TileCache tool for tile caching, i can use it with Leaflet. Like this;

var layername = L.tileLayer.wms("https://tile-cache-server-adress/tilecache/tilecache.py", {
    layers: 'LayerName',
    format: 'image/png',
    transparent: true,
    attribution: "",
    maxZoom: 20
});

So, is there another way to implement tile caching a vector layer (WMS) and publish it as a WMS again for using with Leaflet?


texts: Modifying Speed Profile and daily profile table for historical traffic in ArcMap Network Analyst

texts:
I am using HERE Traffic Pattern and Navstreet data to create a Historical traffic network for my analysis, but as mentioned in this guide: https://desktop.arcgis.com/en/arcmap/latest/extensions/network-analyst/traffic-historical-10-1-and-later.htm
under "How Network Analyst models historical traffic" my daily profile data is not in scale factor (i.e. ratio of average speed at a given point of time and freeflow speed) and hence I am getting a total green (freeflow) traffic map which is incorrect.
I have derived Freeflow Speed (from Navstreet Feature class) and Travel Speed (HERE Traffic Pattern) for given links.
In excel, I used vlookup to fill the data in a new sheet (168 columns as have 1*24 Hour slices for 7 days), divided the average speed by freeflow speed, created a pivot of these exported that into a new table and used vlookup again to fill the new pattern ID. But this process is valid only for 1048576 Road links because of excel limit whereas I have 3360000+ road links.
Now my question is is there any other way to convert my daily profile to scale factors or do I need to create a new table? If so How can I do that without using excel?
I Have attached Screenshot of my Street Profile, Daily Profile and Spatial Join of navstreet-Street profile for better understanding
Thanks!
 Street Profile (Traffic Pattern)
 Daily Profile Table
StreetProfile_Navstreet Spatial Joined (as used in network fgdb)


texts: Extracting donut holes as separate layer in QGIS?

texts:
Does anyone know how I can extract the holes in a polygon layer as a separate layer in QGIS?

An example could be a country with a number of lakes cut out of the polygon - I need the lakes.


texts: How to calculate proportion of vegetation (Pv) for the Landsat 8 image collection?

texts:
I have calculated the Brightness temperature & NDVI for the filtered Landsat 8 collection by making functions.
But I am not able to calculate Proportion of vegetation(Pv) for it. As to calculate Pv we need max & min values of NDVI, so I'm also not able to calculate it. I think that If I'll be able to calculate max. & min values of NDVI then it would be easy to calculate Pv.
Also, if there is any other way to calculate Pv ; i.e. by just putting -1 & 1 in place of minNDVI & maxNDVI respectively in the method to calculate Pv or something like that then also suggest me.
Here is the link to my code:https://code.earthengine.google.com/a0a3671455d9ce767ca3ef34a4a8f9cd


texts: Overlapping map tiles with mapnik's generate_tiles.py script

texts:
I'm trying to generate map tiles using the mapnik's generate_tiles.py script, but the resulting tiles are overlapped.

Example

The tiles have the following directory structure: zoom_level/x_coordinate/y_coordinate.png.

Tiles 3/4/1.png , 3/4/2.png and 3/4/3.png have some countries shown twice (for example: Turkey and a part of Italy) like in the following image:

Begin Edit

Here is how the tiles should look at the same zoom level as the previous image (zoom level: 3; these tiles were generated using MOBAC with the "OpenStreetMap Mapnik" Map source):

As you can see, in this case there aren't any tiles with overlapping content (or duplicate, if you prefer). Another strange thing is that only the tiles from the middle of the map (3.png and 4.png) look similar (cover approximately the same area). The other tiles (0.png , 1.png and 2.png) look totally different and I don't understand why!!!

Note: this is the first time I generate my own tiles so it's very likely that I'm doing a stupid mistake. It might be related to the projection I use (see below for the osm.xml).

End Edit

I tried setting the self.m.buffer_size to 0 (zero) (like this page says; anyway, that's regarding labels, which I don't need) in:

def render_tile(self, tile_uri, x, y, z):
    # ...

    render_size = 256
    self.m.resize(render_size, render_size)
    self.m.zoom_to_box(bbox)
    self.m.buffer_size = 128

    # Render image with default Agg renderer
    im = mapnik.Image(render_size, render_size)
    mapnik.render(self.m, im)
    im.save(tile_uri, 'png256')

Could it be a problem with the projection? Here are the important bits from my osm.xml file:

<Map bgcolor="#f2efe9" srs="+proj=latlong +datum=WGS84">

and

<Layer name="countries" status="on" srs="+proj=latlong +datum=WGS84">

Can someone help me? I already spent almost a week searching and documenting myself about this process and I got stuck here.


texts: Shapefile not showing up in QGIS, and doesn't give coordinates in (Lat, Lon) in R

texts:
I am attempting to find the coordinates associated with the Traffic Analysis Zones from this survey.

When I open it in QGIS it shows up on the sidebar as a point, not a block, and when I open the shapefile in R it gives coordinates like this:

coordinates(shape2)

       [,1]    [,2]
0    2406734 1631930

1    2433452 1626985

2    2447257 1617579

3    2391080 1621358

4    2420945 1617582

5    2373529 1616181

6    2410228 1612591

Here is a link: http://www.surveyarchive.org/

Column: Raw Data, Row: Atlanta, then click on 2011

Also, I've realized it doesn't come with a .prj file. Is there a fix for this?


texts: Trying to use the GeoWebCache "gmaps" service

texts:
I am using GeoServer 2.4.5 with built-in GeoWebCache 1.5.0.

I'm integrating GeoServer with a Google Maps V3 api frontend, and I got it to work reasonably well using the WMS service.  However, I noticed that there's a builtin 'gmaps' service, as described here: http://geowebcache.org/docs/1.5.0/services/gmaps.html

I've tried to use it, with a url like:

http://myserver:8080/geoserver/ok/gwc/service/gmaps?layers=ok:OK_Twp&zoom=9&x=114&y=200&format=image/png

however, I get the following response:

This XML file does not appear to have any style information associated with it. The document tree is shown below.
<ows:ExceptionReport xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0"     xsi:schemaLocation="http://www.opengis.net/ows http://myserver:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="InvalidParameterValue" locator="service">
<ows:ExceptionText>No service: ( gmaps )</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

For some reason, it can't find the 'gmaps' service.  Any ideas?  I've made sure to enable "Global Services", but can't find any configuration regarding the gmap service in the GeoServer web interfact.


texts: Cannot open MapInfo created .TAB files in QGIS

texts:
Newly installed QGIS and trying to open a .TAB file created with MapInfo Pro but fails to create shapefile receiving message ".TAB file not supported".

I Followed instructions > layers > add vector layer >.

Any suggestions please?


texts:
I have a Tile server developed in JavaScript using express. Here you can see the code.

But I would like to know if its possible to load a cluster layer created using Mapbox and tiles to cluster layer. With GeoJSON file works perfectlty but with tiles, render all circles but circle-color and circle-radius doesnt work.

All circles are equal and black by default.

My code:

<script>
        mapboxgl.accessToken = '<your access token here>';
        var map = new mapboxgl.Map({
            container: 'map',
            style: 'mapbox://styles/mapbox/dark-v10',
            center: [-103.59179687498357, 40.66995747013945],
            zoom: 3
        });

        map.on('load', function () {

            map.addSource('clustersource', {
                type: 'vector',
                tiles: ['http://192.168.240.71:8090/tile/terremotos_3e7010c6_2ac0_4895_a392_8dfe212c113f/{z}/{x}/{y}.pbf'],
                cluster: true,
                clusterMaxZoom: 14,
                clusterRadius: 50
            });

            map.addLayer({
                id: 'clusters',
                type: 'circle',
                source: 'clustersource',
                'source-layer': 'terremotos_3e7010c6_2ac0_4895_a392_8dfe212c113f',
                paint: {
                    'circle-color': [
                        "step",
                        ["get", "point_count"],
                        "#51bbd6",
                        100,
                        "#f1f075",
                    ],
                    'circle-opacity': 1,
                    'circle-radius': [
                        'step',
                        ["get", "point_count"],
                        20,
                        100,
                        30
                    ]
                },
                layout: {
                    visibility: 'visible'
                }
            });

        });
    </script>


texts: Paths layers disappear on map after saving in QGIS 3

texts:
I am quite new to GIS. I am working on a historical map with three sets of points. When I create a path between those points with the "Points to Path" tool (via Ctrl+Alt+T), it shows up nicely.

However, when I save the map and open it again later, the paths layers are still visible in the layer toolbar, but no longer visible on the map, not even when I click "Zoom to map layer". I am working in QGIS 3.2 (Bonn).

How can I fix this?


texts: How could we get presentation of FOSS4Geurope 2017?

texts:
I was to the FOSS4Geu 2017 and I could not see all the presentation I would to see. So, as the great team try to get all of them, I thought it's possible to get it throught the website.
Unfortunatly I can't figure out where it's possible to download them.
Do you know where presentations are available ? maybe  on a specific FOSS4G page dedicate to that (which I never earded talking about before)!


texts: Can't move some vertices when editing in QGIS 3

texts: I have a shapefile in QGIS and I can move some of the vertices using the vertex tool but when I try the same approach on another vertex, it does not work.


texts: Return single image from image service (mosaic dataset)

texts:
I have managed to publish multiple raster (tiff) under one mosaic dataset as Image Service.

I also manage to use ArcGISImageServiceLayer to load the ImageServer above.

but how do I retrieve a single raster image based on their name?

I tried to set the ImageServiceLayer initialExtent and fullExtent, but ImageServiceLayer always use the Map Extent which always return all the rasters.

I also used "setDefinitionExpression" method in ImageServiceLayer but there is nothing being drawn on the map.

The closest that I manage to do was:
1. use query on ImageServer with Name = ''
2. Get the Extent of the return record
3. Create ImageServiceParameters with Extent in #2 above as value
4. use exportMapImage and pass in ImageServiceParameters in #3 above.

When I check Chrome debugging console, #4 above return the correct image that I want to load.

How do i load the MapImage returned by the method callback into current ImageServiceLayer?


texts: Why do zones differ for WGS84 UTM in Google Earth Pro vs. in QGIS 3.4?

texts:
For example, zone UTM zone 12S in Google Earth Pro contains southern UT and AZ...

while in QGIS, zone 12S appears to be in the southern hemisphere. Probably a ridiculous question, but not quite sure what am I missing.


texts: Is it possible to use pgsql2shp utility in PL/pgsql environment?

texts:
I intend to use PL/pgSQL to generate calculation points on streets with respect to address locations. I am aware that using pgsql2shp utility, a shapefile can be generate using postgis query.

I was wondering whether is it possible to use pgsql2shp utility in PL/pgSQL functions?

What I need to do is to:

1) first generate the shapefile of calculation points using shp2pgsql utility on the go in plpgsql

2) then use this shapefile for some geo-processing in the same plpgsql environment for example using ST_DWithin() or ST_Distance() etc.


texts: QGIS default projection used to render a WGS84 map

texts:
This question seems naive at first glance, but it's subtle. I was recently pondering... How does QGIS actually represent maps and features when we select the 'WGS84' CRS at the project level?
I surprisingly cannot find the information here; https://docs.qgis.org/3.16/en/docs/user_manual/working_with_projections/working_with_projections.html (neither in the software itself).
In other words, what is QGIS default 'under-the-hood' projection used when WGS84 is selected for a project?

Fig. 1: WGS84 selected as the project CRS to render an OpenStreetMap background map - my screen is definitely not ellipsoidally shaped.
This question can obviously extend to any other GIS piece of software, so instead of asking for each, if in addition you can tell if it exists a standard which is (almost) always used when displaying WGS84 material on a flat screen, that would be even better to know. It seems that it is not 'embedded' in the WGS84 CRS itself (which is a good thing) as per; https://epsg.org/crs_4326/WGS-84.html?sessionkey=3gcvw501ip and:
$ projinfo epsg:4326
PROJ.4 string:
+proj=longlat +datum=WGS84 +no_defs +type=crs

WKT2_2018 string:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

This still remains a Geographic space for me, but which is in 2D only, so how does it actually goes to my flat screen 2D Cartesian space?
Note; I intuitively think it's using an equidistant cylindrical (Plate Carrée) (a.k.a. equirectangular projection) projection, but I definitely wish to strike out that 'intuitively' word from the current sentence.
More on this projection: https://proj.org/operations/projections/eqc.html
This is a parenthesis, but I also think it would be a nice feature to see in a user-friendly way whether a CRS is a Geographic one or a Projected one. For the moment one has to watch to the WKT string to figure it out:

In my researches, I found this useful article on the ESRI blog:

Source: https://www.esri.com/arcgis-blog/products/arcgis-pro/mapping/gcs_vs_pcs/


texts: QGIS different scale levels in rule-based styling - in one rule

texts:
I have a layer with 4 roadclasses. For this layer I created a rule based symbolization, with 4 rules. In every rule I filter another roadclass and set a different symbolization and scale for this rule.
Now I want to create a fifth rule "road direction", which should show the road direction as an arrow. But this arrow should be visible in different scale levels, depending on the roadclass.

For example:
the roadclass 'A' the arrow should be visible from a scale 1:1500.000
the roadclass 'B' the arrow should be visible from a scale 1:400.000
and so on

Can I set this in one rule? The aim is, that I can activate or deactivate the rule "road direction" with one click. I don´t want to create 4 different rules and the arrow should not be added to the 4 rules for the roadclasses.


texts: Can ST_BUFFER and ST_WITHIN be used for linestrings within a linestring?

texts:
I am working with MariaDB and trying to get data based on the condition of linestrings that are within another linestring. To allow linestrings not having to be exactly within (due to GPX variability in coordinates, see note), I want to include a buffer.
SELECT * FROM tableX WHERE ST_WITHIN(ST_BUFFER(linestring,r),
  (SELECT linestring FROM tableY WHERE id=@id) )

I expect to retrieve at least one row, since this row is also retrieved with the same query without the buffer (so, I would expect this is equal to r=0). However, no matter how small or big r is, I get an empty table.
In the documentation it says that ST_BUFFER applies to all geometry objects, but I have only seen examples with points in either linestrings / polygons. It made me wonder if it is even possible to apply this function to linestrings within a linestring?
note: eventually I want to see if a GPX (blue) covered a specific path (purple).
Due to GPX variability, the path is not acknowledged by ST_WITHIN to be within the GPX track:

Hence, the buffer function would be very useful here.


texts: Is there an easy way to read output from Open Data Kit in QGIS?

texts:
I use Open Data Kit (ODK) to do field collection of points of interest. ODK outputs each collected point as an XML-file with a structure like so:

<?xml version='1.0' ?>
<data id="build_My_POI_Type_1234">
<meta>
<instanceID>uuid:123-abc-456</instanceID>
</meta>
<NAME>Place name</NAME>
<KOORD>56.1234 13.1234 70.1234 7.1234</KOORD>
</data>

Is there a simple solution for reading this as a point feature in QGIS, or do I have to parse it via Python console or something?


texts:
I want to display data (tiles and utfgrid info) inside MSSQL on top of google maps using Mapnik.

I have looked into mapnik and setup a small prototype both with tilestache and also with node-mapnik following http://www.nautilytics.com/blogs/2014/04/15/running-route-tile-server/ .

I prefer to use nodejs and the node-mapnik method, as I think  using nodeJS would be easiest to modify the server-side to allow for dynamic tiles.

I need to use the ogr virtual format to connect to MSSQL in order to get live data to create the tiles, or alternatively get data from a redis server. ( Found NO information about using mapnik and a redis datasource). Connecting to MSSQL is something currently broken in Mapnik, but should be fixed soon:

I am having issues finding a solution that allows for serving different tiles based on a session/SQL query similar to the idea mentioned : http://forumone.com/insights/mapping-county-health-rankings-mapnik-nodejs-and-postgresql/

I have an idea how I can take additional GET parameters in the tileserver such as /{z}/{y}/{z}/png?sql=SELECT id,lng,lat,type FROM LOCATIONS then parse it with nodeJS to extract the live data and serve the points using Mapnik.

I know I can define a datasource on the fly using nodeJS using postgres but not sure if this can be done with OGR.
My hope is that using the Mapnik OGR virtual layer that the query is processed only once, and not every tile request.

I also have absolutely no idea how to get the data to display on the map. I know a good technique is using UTFGrids but I havent been able to find out to translate data inside a database to usable UTFGrids, only tutorials on how to do it using non-dynamic data using TileMill.

So I need help on how to create session-based tiles and how to generate UTFGrids from a database using Mapnik


texts: I have a series of spot NDVI data for Australia and I need to extract some information from it. Unfortunately there is no spatial reference system attached to it. This is one example of the data (HDF files). Do I have the wrong dataset? Is there any other NDVI products useful for Australia?


texts:
Doe is exists a way to remove all the places from the OSM map? Also like shops, bars, restorations, hotels etc.

I wish to use a maps with a less info. Is it possible to do this from the original OSM server? Maybe like a option in URL or something else? I use Leaflet library for my project, maybe some option in it?


texts: SpatiaLite splitting large polygons from imported shape file into multiple rows

texts:
I have this shapefile:
shapes.shp
in which I have imported in SpatiaLite using this command:
spatialite_tool --import --shapefile shapes --db-path db.sqlite --table shapes --charset UTF-8 --srid 27700 --coerce-2d -k

which provides me with a geometry column which seems valid:
SELECT Count(*), GeometryType("geometry"), Srid("geometry"), CoordDimension("geometry")
FROM "shapes"
GROUP BY 2, 3, 4

and the output:
Count(*)    GeometryType("geometry")  Srid("geometry")  CoordDimension("geometry")
----------  ------------------------  ----------------  --------------------------
509282      MULTIPOLYGON              27700             XY

I then transform the geometry to SRID 4326 using the following:
UPDATE shapes SET multi_geom = TRANSFORM(geometry, 4326);

I then run the follwing:
SELECT Count(*), GeometryType("multi_geom"), Srid("multi_geom"), CoordDimension("multi_geom") from shapes GROUP by  2, 3, 4;

Which proves me with this:
Count(*)    GeometryType("multi_geom")  Srid("multi_geom")  CoordDimension("multi_geom")
----------  --------------------------  ------------------  ----------------------------
509282      MULTIPOLYGON                4326                XY

All seems good.
Then I apply the Spatial Index via:
SELECT CreateSpatialIndex("shapes", "multi_geom");

So now I can query quite quickly any items that intersects a bounding box and it works well. The problem I now have is that I would like to display these golygons on a map using AsGeoJSON(), and though it works fine in most cases I do however have some very large polygons:
SELECT ST_NPoints(multi_geom) AS pointCount FROM shapes ORDER BY pointCount DESC limit 10;
pointCount
----------
460061
303778
248892
241198
238410
238153
205812
193641
192667
178179

Ultimately what I would like it to be able to use a form of "ST_Subdivide" function or something similar which will break up the polygons into different "rows" so that I am not sending 20MB of GeoJSON to a browser when a large polygon intersects with my bounding box. I can deal with some overlap, however some of the polygons span 10's of kilometres which in my eyes is wasted bandwidth and compute. I have been trying various different "methods" to try and achieve this to no avail, so I am certain I am missing something fundamental here.
An example of the sheer size of polygon could be seen running this query:
SELECT multi_geom
FROM shapes
WHERE ROWID IN (
        SELECT ROWID
        FROM SpatialIndex
        WHERE f_table_name = 'shapes'
            AND f_geometry_column = 'multi_geom'
            AND search_frame = BuildMbr(
                -2.022857666015625,
                52.109351808806345,
                -1.9955635070800781,
                52.11989325575632,
                4326
            )
            AND ST_Intersects(
        multi_geom,
                BuildMbr(
                    -2.022857666015625,
                    52.109351808806345,
                    -1.9955635070800781,
                    52.11989325575632,
                    4326
                )
            )
);

While in reality I would like the polygon broken up so I don't have to render it all.


texts: Calculating coverage of vector polygon over other polygons?

texts: I have a shapefile file that contains a LiDAR footprint and I'm trying to figure out how much a county is covered by the LiDAR footprint. Two separate shapefiles. The attribute table for the county layer is separated by counties with an area. The LiDAR footprint is one large footprint that covers multiple counties.


texts: How to recursively loop through parent polygons intersections to get smallest (child) polygons with no overlaps?

texts:
I'm struggling with a problem for a couple of days and realized many people also get stuck when the topic is intersections in PostGIS (v2.5).
That's why I decided to ask for a more detailed and generic, common question.

I have the following table:

DROP TABLE IF EXISTS tbl_foo;
CREATE TABLE tbl_foo (
    id bigint NOT NULL,
    geom public.geometry(MultiPolygon, 4326),
    att_category character varying(15),
    att_value integer
);
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
    (1, ST_SetSRID('MULTIPOLYGON (((0 6, 0 12, 8 9, 0 6)))'::geometry,4326) , 'cat1', 2 );
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
    (2, ST_SetSRID('MULTIPOLYGON (((5 0, 5 12, 9 12, 9 0, 5 0)))'::geometry,4326), 'cat1', 1 );
INSERT INTO tbl_foo (id, geom, att_category, att_value) VALUES
    (3, ST_SetSRID('MULTIPOLYGON (((4 4, 3 8, 4 12, 7 14,10 12, 11 8, 10 4, 4 4)))'::geometry,4326) , 'cat2', 5 );

It looks like this:

I want to get all the child polygons based on the intersection of the parent polygons. For the outcome, it'd be expected:

The child polygons with no overlap between them.
An column containing the sum of the value of their parent polygons,
An column containing the count of parent polygons of one category
An column containing the count of another category
An column containing the category of the child polygon, based on the following rule: -If ALL the parent polygons are from one class, the child polygon also have this class. Else, the category of the child polygon is a third category.

So it'd looks like it:

So, in the end, the output table generated (for this example) will have 7 rows (all the 7, non-overlapping, child polygons), containing columns of category, sum_value, ct_overlap_cat1, ct_overlap_cat2

The following code I started, gives me the individual intersections, comparing one parent with another.

SELECT
(ST_Dump(
    ST_SymDifference(a.geom, b.geom)
)).geom
FROM tbl_foo a, tbl_foo b
WHERE a.ID < b.ID AND ST_INTERSECTS(a.geom, b.geom)
UNION ALL
SELECT
ST_Intersection(a.geom, b.geom) as geom
FROM tbl_foo a, tbl_foo b
WHERE a.ID < b.ID AND ST_INTERSECTS(a.geom, b.geom);

How do I recursively loop through the result of this mentioned code, that, independent of the number of overlap polygons I always get its 'smallest' (child) polygons (Fig. 2)?


texts: QGIS freezes after making and saving many edits to a layer

texts:
I work with QGIS 2.4 on windows 7.

The problem occurs when I work on a layer and save it many times. After a random number of saves, pressing the save button once more causes QGIS to freeze, all processor being used for more than two minutes. Exiting the edit mode for the layer also causes that behaviour.

The only solution to the problem is relaunching qgis. Saving the project does not cause a freeze, but clearing the  actual project does. And with that button the freeze will not recover, so I have to kill QGIS.

Has anyone encountered the same issue, or has any idea about it ?


texts: creating Polygon from LineString and a width

texts:
I have to create a Polygon from a lineString and a width. Ex: Given a street as a lineString (Or a list of points) and a width, create the polygon (the road). I am using Geotools, but I couldn't find anything helpful in the API.

The math for doing this manually is requiring to much work. Is there a simpler solution for this kind of problem?

Thanks.


texts: Overlaying lines and aggregating their values for overlapping segments using R?

texts:
I have a number of routes (5 in this
example data)
stored in a
SpatialLinesDataFrame called r. These have different origins but all
converge on the same destination
(Manchester Airport):

Each contains the variable All, representing
data
on the number of people regularly using each route.

How can I merge the lines such that the sections where they overlap are allocated a value of
sum(All) for all the lines passing through there?

The red and yellow lines, for example, have All values of 210 and 395,
so the orange line going north-south should be a segment with a value of
605:

> r_overlap$All + r[4,]$All
[1] 605

But I need all the merged lines in one layer, with each new Lines object representing a segment with a single sum(All) value, not just 1 as above.

The data above is stored at raw.githubusercontent.com/npct/pct-data/master/test-data/airport.geojson  so the solution can be demonstrable and reproducible.
An R-based solution would be preferable to me.

The below code was used to load and view the data using the new geojsonio and leaflet packages:

pkgs <- c("geojsonio", "leaflet", "sp")
lapply(pkgs, library, character.only = T)
download.file("https://raw.githubusercontent.com/npct/pct-data/master/test-data/airport.geojson", destfile = "l.geojson", method = "wget")
r <- geojsonio::geojson_read("l.geojson")

plot(r) # check the data is there
r@data # look at the data

r_overlap <- gIntersection(r[1,], r[4,])
plot(r)
plot(r_overlap, col = "red", add = T, lwd = 5)

r_overlap$All + r[4,]$All

leaflet() %>%
  addTiles() %>%
  addPolylines(data = r, color = c("red", "blue", "green", "yellow", "orange"))


texts: Identifying double/overlapping polygons using ArcGIS Desktop?

texts:
I have two datasets in which the information about the area are deviating:

an observation room with three corridors (A = 25000 ha, B = 30000 ha and C = 20000 ha)
land use dataset provided by an external institution (A = 35000 ha, B= 45000 ha and C = 37500 ha)

I do not know how the land use data set was approached, but I assume that the result of it includes double or overlapping polygons. Can anyone show me a way, how to identify and delete those features t receive the same area as in the observation room.

I use an ArcGIS 10.3 basic license.


texts: Python syntax error -can't figure out why?

texts:
This should be simple. I hope I don't get down votes!

I just got the following code from an article about using Google Maps API and ArcGIS to calculate routes without network analysis and just can't get it to run. I am getting a syntax error on line 7. I have tried to type in other values before it and they are fine, also tried the to change the value "fromFile" to something else and it still gives the same error.

Any ideas?

#Estimating O–D travel time matrix by Google Maps API: implementation, advantages, and implications
#http://dx.doi.org/10.1080/19475683.2011.625977
import arcpy
import urllib
import time
from xml.etree.ElementTree import XML, fromstring, tostring
fromFile = arcpy.GetParameterAsText(0)
toFile = arcpy.GetParameterAsText(1)
Resultfile = arcpy.GetParameterAsText(2)
Result = open(Resultfile,”w”)
Result.write(“FromFID,toFID,TravelTime”)
Result.write(“\n”)
fromCursor = arcpy.SearchCursor(fromFile)
toCursor = arcpy.SearchCursor(toFile)
fromRow = fromCursor.reset()
fromRow = fromCursor.next()
while (fromRow!=None):
fromX = fromRow.shape.centroid.X
fromY = fromRow.shape.centroid.Y
fromFID = fromRow.FID
arcpy.AddMessage(str(fromFID))
toCursor = arcpy.SearchCursor(toFile)
toRow = toCursor.reset()
toRow = toCursor.next()
while (toRow!=None):
toX = toRow.shape.centroid.X
toY = toRow.shape.centroid.Y
toFID = toRow.FID
arcpy.AddMessage(str(toFID))
googletext = “http://maps.googleapis.com/maps/api/directions/xml?origin=(“ + str(fromY) + “,” + str(fromX) + “) &destination=(“ + str(toY) + “,” + str(toX) + “) &sensor=false”
time.sleep(3)
xmlfile = urllib.urlopen(googletext)
xml = xmlfile.read()
value = “NA”
dom = fromstring(xml)
nodelist = dom.getchildren()
if (nodelist[0].text =  = “OK”):
arcpy.AddMessage(nodelist[0].text)
route=nodelist[1]
leg=route.getchildren()[1]
duration = leg.find(“duration”)
value = duration.getchildren()[0].text
else:
arcpy.AddError(nodelist[0].text)
Result.write(str(fromFID))
Result.write(“,”)
Result.write(str(toFID))
Result.write(“,”)
Result.write(value)
Result.write(“|n”)
toRow = toCursor.next()
fromRow = fromCursor.next()
Result.close()
del fromCursor
del toCursor


texts:
Short version

How can I remove SRID information from PostGIS geometry? Is ST_SetSRID(geom, 0) well-defined?

Long version

I am rendering maps, using PostGIS. To do this, I create a temporary table, and populate it by ST_Transform():ing my geometries from EPSG:4326 to my desired map projection (usually EPSG:3006, SWEREF99TM) and using ST_Affine() to scale and translate them into map coordinates (90-per-inch units on the paper). I then do some processing on this temporary table (cartographic generalization of different kinds), before marshalling it all out to the styling stage.

Apparently, the result of ST_Affine(ST_Transform(..., srid), ...) is considered to be “in” that SRID, even though the ST_Affine makes this dubious, at best. This means that whenever I want to add geometry from other sources into the mix, that, too, has to be tagged with that SRID, even though it is actually in “paper space”. This means that I have to keep track of the SRID in various places in my program, which are really only concerned with things like “25 mm from the left side of the paper”. If I don’t, PostGIS complains that I am doing
operations on geometries with different SRID:s.

It seems like PostGIS happily performs things like ST_Intersects on geometries with no SRID, so if I could get everything to this state, I would be fine.


texts: Plugin that selects closest point features from a mouse click

texts:
I'm trying to build a plugin that would select let's say 8 closest points from a mouse click. I've found a working plugin that selects the closest feature from a mouse click here: https://www.lutraconsulting.co.uk/blog/2014/10/17/getting-started-writing-qgis-python-plugins/

So far I've managed to make it select only point features and I know how to make it select 2nd or 3rd closest point, but I have trouble making it select more features. The following doesn't return any selection:

class NearestFeatureMapTool(QgsMapTool):

    def __init__(self, canvas):

        super(QgsMapTool, self).__init__(canvas)
        self.canvas = canvas
        self.cursor = QCursor(Qt.CrossCursor)

    def activate(self):
        self.canvas.setCursor(self.cursor)

    def canvasReleaseEvent(self, mouseEvent):
        """
        Each time the mouse is clicked on the map canvas, perform
        the following tasks:
            Loop through all visible vector layers and for each:
                Ensure no features are selected
                Determine the distance of the closes feature in the layer to the mouse click
                Keep track of the layer id and id of the closest feature
                Select the id of the closes feature """

        layerData = []

        for layer in self.canvas.layers():

            if layer.geometryType() != QGis.Point:

                # Ignore this layer as it's not a vector
                continue

            if layer.featureCount() == 0:
                # There are no features - skip
                continue

            #layer.removeSelection()

            # Determine the location of the click in real-world coords
            layerPoint = self.toLayerCoordinates( layer, mouseEvent.pos() )

            shortestDistance = float("inf")
            closestFeatureId = -1

            # Loop through all features in the layer
            for f in layer.getFeatures():
                dist = f.geometry().distance( QgsGeometry.fromPoint( layerPoint) )

                shortestDistance = dist
                closestFeatureId = f.id()

                info = (layer, closestFeatureId, shortestDistance)
                layerData.append(info)

        if not len(layerData) > 0:
                # Looks like no vector layers were found - do nothing
            return

            # Sort the layer information by shortest distance
        layerData.sort( key=lambda element: element[2] )

        selected_fid = []
        for closestFeatureId in layerData:
            if len(selected_fid) != 8:
                selected_fid.append(closestFeatureId)
            else:
                break
        layer.select(selected_fid)


texts: Reprojecting EPSG:3412 to EPSG:4326 using GDAL?

texts:
I want to reproject "South polar stereographic" projection to "gcs WGS 84". I have tried the solution given to Re-projecting from Polar projection to EPSG:4326?

through this commands I am getting output in format of WGS84. Moreover latitude value is also correct but longitude values are surprisingly different from expected. what are all possible reasons for this? Code is provided below

gdal_translate -a_srs "+proj=stere +lat_0=-90 +lat_ts=-60 +lon_0=100 +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs" -a_ullr 3818607 -4437313  -4944393  4579687 -of VRT input.tif temp'.vrt
gdalwarp -overwrite -t_srs EPSG:4326 -te -180 -90 180 90 temp.vrt projected.tif


texts: Converting image to shapefile using Google Earth Engine

texts:
I have performed a supervisioned classification to identify lakes in a study region.
How can I transform a classification that are in Image to ShapeFile?
The code that I wrote:
var image = ee.Algorithms.Landsat.simpleComposite({
  collection: l8.filterDate('2015-01-01', '2019-12-31'),
  asFloat: true
});

var image = image.clip(geometry);

// Use these bands for prediction.
var bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B10', 'B11'];

// Manually created polygons.

// Make a FeatureCollection from the hand-made geometries.
var polygons = ee.FeatureCollection([
  ee.Feature(nonlake, {'class': 0}),
  ee.Feature(nonlake1, {'class': 0}),
  ee.Feature(lake1, {'class': 1}),
  ee.Feature(lake2, {'class': 1}),
  ee.Feature(lake3, {'class': 1})
]);

// Get the values for all pixels in each polygon in the training.
var training = image.sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: polygons,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 30
});

// Create an SVM classifier with custom parameters.
var classifier = ee.Classifier.libsvm({
  kernelType: 'RBF',
  gamma: 0.5,
  cost: 10
});

// Train the classifier.
var trained = classifier.train(training, 'class', bands);

// Classify the image.
var classified = image.classify(trained);

//Map.addLayer(polygons, {}, 'training polygons');
Map.addLayer(classified, {min: 0, max: 1, palette: ['white', 'black']},
             'delakeation');


texts: Why do temporary features disappear when hovered over?

texts:
With DrawFeature I'm putting new features to my map. With the hover option of SelectFeature it is possible to get feature information when hovering over it.
Now I have the problem of disappearing features on the map when hovering over newly added (=temporary) features?

Any solutions how to style the features so that they stay on the map during storage/editing process?

Here is the relevant code I use styling the point features:

var defaultPoint = new OpenLayers.Style({ label: "D", pointRadius: 8, fillColor: "red", fillOpacity: 0.7, strokeColor: "black", strokeWidth: 3});
var selectPoint  = new OpenLayers.Style({ label: "S",pointRadius: 8, fillColor: "blue", fillOpacity: 0.7, strokeColor: "black", strokeWidth: 3 });
var tempPoint    = new OpenLayers.Style({ label: "T",pointRadius: 8, fillColor: "green",    fillOpacity: 0.7, strokeColor: "black", strokeWidth: 3 });

var stylemapPoints = new OpenLayers.StyleMap({'default': defaultPoint, 'select' : selectPoint , 'temporary': tempPoint});

This is the code for the Hover control:

HoverSelectControl = new OpenLayers.Control.SelectFeature([layerName],
                        { hover:true,
                          highlightOnly:true,
                          renderIntent: "select"
                         });
HoverSelectControl.handlers.feature.stopDown = false;
map.addControl(HoverSelectControl);
HoverSelectControl.activate();


texts:
I have installed the new Geoserver 2.8.0 and to avoid CORS related errors I wanted to implement a reverse proxy, instead of setting response headers.

In version 2.1.3 there was a <filter-name>Reverse Proxy Filter</filter-name> that you could enable, in the web.xml file.

In version 2.8.0, I could not found such filter. How can I implement a reverse proxy, so I can get rid of CORS related issues?


texts: From several polygons of one layer to several layers/files

texts:
I'm new in QGIS and I have a question. I am sampling an onion crop from a geottif:

I need to export each single rectangle/polygon into a single tiff file. Is there a plug-in for doing so?


texts: QGIS 2.0 OpenLayers Google Maps/Satellite do not work

texts:
I have been using QGIS version 2.0 with Google Satellite layers for a few months. The last time I used was in mid Oct 2014. In mid Nov 2014, Google Satellite or Google Map suddenly does not work any longer. I'm unable to zoom in or zoom to the layer extent. It is showing a static graphic of Google map in a high scale. I tried uninstalling and re-installing the OpenLayers plugin, but it still didn't work. I tried to change the CRS Layers from EPSG:3857 to many different layers. It still does not work.  Can someone please help?

Below is snapshot of what happened to the Google Satellite or Map.


texts: GRASS Toolbar for QGIS missing the tools I expect to see

texts:
I am trying to learn GRASS and QGIS using the instructions here: https://docs.qgis.org/2.6/en/docs/user_manual/grass_integration/grass_integration.html. I am running on Windows and open the version of QGIS with GRASS. When I enable the GRASS plugin, however, I see only two buttons on the toolbar. The panel is telling me to open a GRASS mapset but the button to create a new mapset is not visible in the toolbar. What am I missing here?

I'm using QGIS 2.14.3 with GRASS 7.0.4


texts: Assign raster cells presence/absence values depending on whether polygon overlaps with cells

texts:
I have a raster of continuous numeric values. A certain portion of said raster overlaps with a polygon of interest. I want to know which cells overlap, and assign values to the raster's attribute table corresponding to whether or not there is overlap (e.g. 0 for no overlap with polygon, 1 for overlap).

Ideally it would be great to be able to set an overlap threshold so that only raster cells where there is more than 50% overlap are assigned a "1". Not a must, but would be very helpful.


texts:
I'm trying to run Mort Canty's http://mcanty.homepage.t-online.de/ Python iMAD implementation on bitemporal RapidEye Multispectral images. Which basically calculates the canonical correlation for the two images and then substracts them. The problem I'm having is that the images are of 5000 x 5000 x 5 (bands) pixels. If I try to run this on the whole image I get a memory error.

Would the use of something like pyTables help me with this?

What Mort Canty's code tries to do is that it loads the images using gdal and then stores them in an 10 x 25,000,000 array.

# initial weights
wt = ones(cols*rows)
# data array (transposed so observations are columns)
dm = zeros((2*bands,cols*rows))
k = 0
for b in pos:
band1 = inDataset1.GetRasterBand(b+1)
band1 = band1.ReadAsArray(x0,y0,cols,rows).astype(float)
dm[k,:] = ravel(band1)
band2 = inDataset2.GetRasterBand(b+1)
band2 = band2.ReadAsArray(x0,y0,cols,rows).astype(float)
dm[bands+k,:] = ravel(band2)
k += 1

Even just creating a 10 x 25,000,000 numpy array of floats throws a memory error. Anyone have a good idea of how to get around this?


texts: can I get map's bounding box at any time, after user clicks button

texts:
I use openlayers. I am creating a form that the user can insert a bounding box to the database.

In my page I have the form and a map, so the user can preview what he entered.

So I thought it would be nice if there was also the other way around :

User hits a button like "Get current bounding box", and that gets the bounding box from the map at any given time and transport it to the form. The numbers of the bounding box that is.
I'm talking about the bounding box that the users see, not the whole map's box.

After that, user can hit "Save" and submit the form, or can clear the form, or edit it, or zoom/pan the map and hit "Get current bounding box" again.

I think you got the idea...

Is that possible?

How can I get map's bounding box at any given time?

Thanks in advance


texts:
I've downloaded some USGS DEMs from the National Map, which appear to be .adf files.  I can open them in ArcMap 10.2, but I need to merge them into a single mosaic.  To do this, I've created a raster catalogue inside a geodatabase file, but when I try to 'Load rasters from workspace', the files are not visible.

Not sure what to do next.  The files are clearly Esri-compatible (I can load and view them individually), but the steps to create a mosaic seem unable to read the files.


texts:
I have a dataset of well pumps, and I am trying to geolocate them. Each of the wells has as "Service Address" field, but instead of giving an address, it lists something that looks like this: "NW NE SW 27 22 23" or "SW SW SW 10 23 23" or even "NW NW NW NW 7 23 21", though maybe that last one is a typo.

Are these spatial coordinates of some type?

It seems a bit like an Aliquot part, but the fact that it's (usually) 3 directions and 3 numbers make me think it must be something different.


texts: How to interpret what each return means from an OSR function?

texts:
I have searched internet for the documentation on what different returns from OSR functions like osr.SpatialReference().IsProjected(), osr.SpatialReference().Validate() means but I couldn't find anything.
I am getting 1 or 0 from IsProjected and 0 or 5 from Validate. I want to know what each integer means.
I have this Python code
srs = osr.SpatialReference()
srs.ImportFromESRI([line])
check = srs.Validate()
print(check) # returns 0 or 5

Is there a good source of descriptive documentation for osr/ogr?


texts:
I currently am using the code below to create a List of existing Domains on a Geodatabase, and if there are none to create one. The code does create the domain, but is not creating the list of Domains if they already exist. Did I set up the if statement incorrectly or do I need another function to give me the list of current domains?

I have referenced the following question, How to check if domain already exists? but the print function is not working for me.

domains=arcpy.da.ListDomains(outputGDB)
if len(domains) == 0:
    domName = 'DomStatus'
    arcpy.CreateDomain_management(r'X:\GIS_PROJECTS\City.gdb', domName, '', 'TEXT', 'CODED')
    arcpy.AddCodedValueToDomain_management(r'X:\GIS_PROJECTS\City.gdb', domName, 1, 'Boundary Used')
    arcpy.AddCodedValueToDomain_management(r'X:\GIS_PROJECTS\NATIONAL\ParkServe\Documentation\Test_HW\ParkScore_Loading\New_Mexico__Albuquerque.gdb', domName, 2, 'Added')
    arcpy.AddCodedValueToDomain_management(r'X:\GIS_PROJECTS\City.gdb', domName, 3, 'Deleted')
    arcpy.AddField_management(newFC, 'Status_2017', 'TEXT',field_length=50, field_domain=domName)
    print "Domain Created"
else:
    print domains


texts: How do I set up the R Script folder in QGIS 2.8.1 Processing Toolbox (Advanced Interface)?

texts:
I am currently lacking the R Script folder in my processing toolbox under the advanced interface tab. I suspect this could be a result of the firewall in my company's wifi.

I already have R installed, just need to figure out how to link it. Preferably a method that does not require administrative rights.


texts: ArcGIS Server REST API Query - Where date strings are greater than

texts:
I am working with an ArcGIS Sever 10.41 rest service that contains many scientific records. These records have a date field called sample_date but it's stored as a string. See below.

sample_date ( type: esriFieldTypeString , alias: sample_date , length: 8000 )

I've been trying to query the last year's records to no avail. Since this field is stored as a string I was trying to just get those records which contained the string 2018 but that did not work and returned empty records. See below.

/arcgis/rest/services/folder/MapServer/2/query?where=sample_date+LIKE+'2018'&outFields=*&f=JSON"

Any suggestions for a better approach to querying all records from the last 12 months even though the field is stored as a string?


texts: How to perform a buffer analysis in PostGIS?

texts: Need to Implement Spatial Analysis in Postgresql. Can Any one guide me to accomplish the spatial analysis eg) 2 km buffer for road network layer.


texts:
I try do load WFS layer to my openlayers application according to this tutorial enter link description here
Features from WFS layer are displaying on the map, but i discover that they are not belong to the Vector Layer, which i've created. This layer is empty. What did i wrong?
// Polygons
function styl_warstwy() {
    return [
        new ol.style.Style({
            stroke: new ol.style.Stroke({
                color: [ 51, 165, 255 ,1],
                width: 2,
                //lineDash: [4,8],
                //lineDashOffset: 6
                }),
            fill: new ol.style.Fill({
                   color: 'rgba(51,165,255,0.25)'
                })
            })
    ];
}

var vectorLoader = function(extent, resolution, projection) {
        var url = 'http://mapy.fotoraporty.pl:8080/geoserver/archiwalne_webGIS/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=archiwalne_webGIS:gmina_result&maxFeatures=2500&outputFormat=application%2Fjson' +
            '&srsname=EPSG:4326&bbox=' + extent.join(',') + ',EPSG:4326';
        $.ajax(url).then(function(response) {
            var format = new ol.format.GeoJSON();
            var features = format.readFeatures(response,
            {featureProjection: projection});
            source.addFeatures(features);
        });
    };

var loadFeatures = function(response) {
        var features = vectorSource.readFeatures(response);
        vectorSource.addFeatures(features);
      };

var vectorSource = new ol.source.Vector({
        loader: vectorLoader,
        strategy: ol.loadingstrategy.bbox
      });

 var vector_kml = new ol.layer.Vector({
    source: vectorSource,
    title: "Gminy",
    style: function(feature) {
            style.getText().setText(feature.get('JPT_NAZWA_'));
            return style;
    },
 });


texts: PostgreSQL cannot work when opens SSL on Windows 10

texts:
I am not sure this topic will fit here but I don't know ask it where.
I am trying to open the SSL of PostgreSQL 10.16 on Windows10.
I read a lot of documents about creating an SSL certificate to PostgreSQL but documents spend on Windows 10 are rare and don't detail.
These are the steps which I did:
Step 1, I download the OpenSSL version of Windows from https://slproweb.com/products/Win32OpenSSL.html
and install it with the path C:\OpenSSL-Win64, setting system variable.
Step 2, I use cmd.exe run as an admin to create server.key with command-line:
  genrsa -out server.key 4096

,  set appropriate permission and owner on the private key file (here https://stackoverflow.com/a/51463654)
  icacls server.key /reset
  icacls server.key /inheritance:r /grant:r "CREATOR OWNER:F"

continue to create the server certificate:
 req -new -x509 -days 1826 -key server.key -out server.crt

Step 3,
Since I am self-signing, I use the server certificate as the trusted root certificate, so I have 3 files: server.key, server.crt, and root.crt (this is a copy of server.crt)
I cut these three files to C:\Program Files\PostgreSQL\10\data
Step 4, I am setting postgresql.conf:
listen_addresses = '*'
port = 5432
ssl=on
ssl_cert_file = 'server.crt'
ssl_key_file = 'server.key'
ssl_ca_file = 'root.crt'

and add command-line to the end of pg_hba.conf:

Finally, I get an error as below when I restart my PostgreSQL

and in my log of PostgreSQL:

I suspect PostgreSQL did not read the 3 files that I put in its data directory.
I have referenced these documents
https://www.howtoforge.com/postgresql-ssl-certificates
https://stackoverflow.com/questions/43959324/chmod-og-rwx-server-key-in-windows
https://www.postgresql.org/docs/10/ssl-tcp.html
I have been messing with it for many days and I don't know how to solve this problem.


texts: Adding point layer in QGIS from SQL Server failing

texts:
I just downloaded QGIS a few days ago (version 2.8.2-Wien) but I am not able to get it to work as I thought it would.  I tried pulling in spatial data from sql server (a "geometry" column containing points).  I click the button that says "Add MSSQL Spatial Layer" and find my point table but when I click the "Add" button it gives me this error (after a minute):

dbname='myDb' host=localhost\sqlserver2014 srid=4326 type=POINT table="dbo"."Locations_Sample"
(spatial_data_point) sql= is an invalid layer and cannot be loaded.  Please check the message log for further info.

The "log message" words in that error were a link in QGIS so I clicked it and it opened a "Log Messages" panel and on that panel (on the "General" tab), it just had the same error message I typed above (with one inconsequential text difference... it says "is an invalid layer - not loaded" instead of "and cannot be loaded").

I tried making that source table ("Locations_Sample") smaller (it was around 120,000 points so I tried just using 100 points) and that smaller table imported with no problems.  I then tried increasing the amount of data in that table and hitting the "Refresh" button in QGIS to find what the row limit was but it succeeded in refreshing every time I tried all the way up to the full 120,000 rows (showing all the points on the map as far as I can tell).  This may seem like a workaround but I don't think it really is because every time I close and re-open the project, it errors when trying to load that layer...also other things (like the heatmap raster generator) are failing possibly due to this "invalid layer" issue.

FYI, I subsequently tried loading the table with 100 records at first and then putting all the records in the source table "Locations_Sample" and that always worked.  Also I think possibly the issue with the heatmap raster generation is that it only works with the initial records I load (100 points).

How can I figure out what the problem is with this layer (and ideally fix it)?


texts: USA Albers Equal Area Conic USGS version (LANDFIRE) projection not read by R?

texts:
I am working with the North America Forest Dynamics dataset. As mentioned it uses the USA Albers Equal Area USGS version projection, which is primarily different from the USA Albers Equal Area projection in the measure of the latitude of origin (former is 23 and latter is 37.5).

The datum is clearly read by Arc as D_NorthAmerica 1983. which is the datum. However, this datum information is not read by R (using raster package). Hence, resulting analysis and written rasters have unknown datum when processing is done in R.

Interestingly, the famous spatial reference website does not have a Proj4 string info for this, which I could manually input into the crs() command. But on the same website, the .Prj file can be downloaded and on opening this ile as a text pad file a datum i.e. North America 1983 datum and hence transformation is mentioned.

Why is R not recognizing the datum around which the transformation has to happen for this projection?