texts: Edit specific indexes from a feature of Attribute Table QGIS

texts:
I have an attribute table with many columns. One of them describes a unique code. The column's name is "specialcode" and the column's type is String(254). I want to do the following:

Sort every feature with exactly 6 indexes in the "specialcode" column and select them.
For every feature from the selection, remove the first 4 indexes, and in the end, only the last 2  indexes (digits) remain as a value for the column "specialcode"

I am assuming it's pretty easy.


texts:
I have a ArcGIS Pro model which takes a polygon shapefile (Site_boundary_master) and spatially joins various layers to it, finally exporting the final output to an Excel spreadsheet:

All the joining layers are kept in various File Geodatabases.

I saved the model as 'AutomationModel' (but in the properties dialogue for it, it's called 'Model222'... wtf?) in a toolbox called 'Automation.tbx'.

Running the model from ModelBuilder works fine.

I then created a .py python file (from a blank .txt file) and wrote a script to run this model this way:

import arcpy
print ('Script started')
# import the toolbox
arcpy.ImportToolbox(r"C:\Automation\Automation.tbx")
print ('Toolbox imported')

#import the model
arcpy.Automation.Model222()
print ('Model imported')
print ('Script finished')

Pasting this code into the ArcGIS Pro python console, it work fine. Again.

However, when I edit the .py file in IDLE (ArcGIS Pro) and then run the script, it doesn't work at all:

Can someone help me fix this?

Is it because I need to tell the script where all the input layers are stored?

If so, how would I write this?

I'm still a novice with python.


texts:
I've tried all of these solutions (How to clip a line layer based on polygons with QGIS?) and still can't clip a shapefile. I'm trying to create a coastal shapefile for one region (NW France) from a world coastline. I created a custom rectangle to use as the clip extent and used vector> geoprocessing> clip to attempt the clip. Still no progress. Has anyone had success using other techniques or can anyone see a problem with my method?

I'm using QGIS 2.12.


texts: Big data rastering using ArcGIS Platform?

texts:
I have the daunting task doing raster analysis on a record set that runs into quite a few millions, and worse still is that the data is in Hive.

So far I have managed to use record set of 2 million and create a raster but, in doing so, I consumed about 15Gb of my server's RAM and took around 5 minutes to complete the whole process. I am expecting the volume of data to increase a few fold.

I have gone through the big data analysis pages of ArcGIS but that did not help me.

Can somebody suggest a better option of doing big data analysis (perhaps by not fetching the data to the code, rather the reverse)?


texts: Is a Nested Loop inevitable in ST_Intersects in PostGIS?

texts:
I'm using PostGIS 9.6 and I have two tables, usgovt and risk. Both have spatial indexes and the database has just had a VACUUM ANALYSE.

usgovt has 142,598 rows, with a mean vertex count of 402 (and a max of 1393536):

category  | integer                     |
geom      | geometry(MultiPolygon,5070) |

risk has just 1 row with a vertex count of 94223:

 geom   | geometry(MultiPolygon,5070) |

I am querying the intersection of the two, but it is slow (380 seconds):

EXPLAIN ANALYZE
  SELECT g.category AS category
  FROM usgovt g, risk r
  WHERE
    ST_Intersects(g.geom, r.geom);

Here is the explain.

 Nested Loop  (cost=4.39..64.15 rows=48 width=2) (actual time=164.231..380784.381 rows=48456 loops=1)
   ->  Seq Scan on risk r  (cost=0.00..1.01 rows=1 width=32) (actual time=0.005..0.006 rows=1 loops=1)
   ->  Bitmap Heap Scan on usgovt g  (cost=4.39..63.09 rows=5 width=34) (actual time=164.223..380718.542 rows=48456 loops=1)
         Recheck Cond: (geom && r.geom)
         Filter: _st_intersects(geom, r.geom)
         Rows Removed by Filter: 71027
         Heap Blocks: exact=23099
         ->  Bitmap Index Scan on usgovt_geom_idx  (cost=0.00..4.39 rows=14 width=0) (actual time=24.500..24.500 rows=119483 loops=1)
               Index Cond: (geom && r.geom)
 Planning time: 0.328 ms
 Execution time: 380855.965 ms

It looks like the slow part is the Nested Loop. Is this inevitable, or is there any way I can avoid it?

(If this seems like an odd query - I want to calculate the area of the intersection, but right now that query doesn't complete. So I'm starting by debuggging whether there's any way to speed up the intersection.)


texts: How to union linestring data distributed on multiple layers?

texts:
I have four shapefiles. Each contains a linestring as shown in the image. Each shapefile is represented as its own layer in QGIS.

I would like to union these linestrings into one polygon which describes the outline of all four shapes as shown in the following image I produced with Gimp. I could not dissolve the parts since they are spread on four layers.

I uploaded the shapefiles for you to try yourself. The archive contains the four separate shapes as well as the merged geometries (Pfinztal-Merged-4326.zip) I produced.


texts: Offline geocoding and reverse geocoding for android

texts:
I'm developing an android application that cannot use internet and need to deal with maps and positioning.

Download OpenStreetMaps data seens to be a good solution to access map information offline. However, I need some basic algorithms running locally on the cellphone, like geocoding and reverse geocoding.

Does anybody know a library for android that can consume offline openstreetmap data and provive geocoding and reverse geocoding on the cellphone?


texts:
I am using the leaflet.pattern library at https://github.com/teastman/Leaflet.pattern

and I want to change the color of my stripes based on a value in the feature just like I do in the following using a getColor Function as part of the style.

function getColor2(d) {
    return d == 9 ? '#d73027':
           d == 8 ? '#f46d43':
           d == 7 ? '#fdae61':
           d == 6 ? '#fee08b':
           d == 5 ? '#ffffbf':
           d == 4 ? '#d9ef8b':
           d == 3 ? '#a6d96a':
           d == 2 ? '#66bd63':
           d == 1 ? '#1a9850':
                '#1a9850';
}

However since the pattern is already defined I am not able to change the color on the fly based on the feature value.

See example at http://jsbin.com/vepiru/1/edit?js,console,output


texts: Conditional styling in FME ChartGenerator

texts:
A vector data set represents different wilderness qualities of district areas (V, 1 and 2) and one human activity zone (N). Our cartography for the areas gives us the following map:

I am using FME ChartGenerator to generate charts for all districts. FME produces excellent charts, with the exception of color scheme selection. I would like to control the pie level color conditional to the wilderness category so that the pies are visually linked to the cartography.

Is there a way to do this which does not involve post-processing color substitution?


texts: Trigger converting multilinestrings to curves

texts:
I created a small web app based on PostGIS / Geoserver / Leaflet.
I can draw polylines from my interface thanks to the Leaflet Draw plugin and send them to PostGIS with a Php and Ajax script.
When my polyline is inserted into PostGIS, a trigger is fired which updates the geometry in my PostGIS table from the received coordinates.
My current problem is that I would like that when I draw a polyline from my interface, my trigger would convert that polyline to a curve as well.
So I looked towards ST_LineToCurve and ST_ForceCurve.
Here is the script in my trigger:
update schema.table SET geom = ST_Transform (ST_GeomFromText ('MULTILINESTRING ((' || coords || '))', 4326), 2154) WHERE geom isnull;

I tried in PostGIS:
SELECT ST_ForceCurve (ST_Transform (ST_GeomFromText ('MULTILINESTRING ((' || coords || '))', 4326), 2154)) FROM schema.table

Which returns geometries in my table, but which displays "unsupported geometries not rendered" in the geometry viewer.
When I do :
SELECT ST_AsText (ST_ForceCurve (ST_Transform (ST_GeomFromText ('MULTILINESTRING ((' || coords || '))', 4326), 2154))) FROM schema.table

the geometries returned in my table are indeed in the form "multicurve".
Now if I add St_ForceCurve in my trigger, nothing is sent to my table.


texts:
In QGIS 3.4 I'm creating a set of ten maps (layouts) from one QGIS project, and each map centres over a particular polygon. In each layout I shift the map to centre it on the page or to fit as best as possible. From that polygon layer I create labels, and I then get one label for each polygon. I can use the tools to manipulate these labels as required, as the labels need to be moved a bit here and there.
  Is it possible to select a label and duplicate that label and have the duplicate locked in the map to another place? This is so I can move the duplicate to the best place for individual map, keeping the original where it was for the previous map. What I can do in ArcGIS: create an annotation layer from labels and then have full control of over that text (copy, resize etc) and have the text locked to the map where I want it - Well I would like to do that in QGIS as easily as I can in Arc. Is this possible?


texts: Calculating area of voronoi polygons from real life data using QGIS?

texts:
I have an area, which I divided into voronoi polygons using data points. I have no problem creating the polygons.

Now the data points are on a coordinate system, where every degree represents one meter. I know try to use the $area command to get the area, but the answer is way to large.

To clarify my problem I will try to give a bit of an example:

This image gives the spatial division of a set of data points, boreholes in this case. They all influence a certain area. This area is defined as their voronoi cell. I have created these cells:

What I now want to know is the area of the polygons. I have tried to use the command $area, but this gives an area which is way too big. I know this because I  did some basic area calculations on some of the simpeler polygons. For example the area for the triangular polygon in the bottom left (BH1101) I calculated to be about 695,000 m2. However using the $area command give me an area of 8.6140950578267e+15 which is in no unit I know equal to what I calculated.

Why is this and how can I fix this?

The data comes from an excel sheet which I converted into a .txt file and put into Qgis. The location is somewhere in Germany, if need be I can get a more precise location, but it shouldnt matter for the results. I run Qgis on the standard settings, using the CRS (WGS84:4326).


texts: Copy query feature layer into a feature class - error 000732

texts:
I am trying to create a python script to copy data from a MSSQL view (includes a geometry field) into a feature class.
They are both held within the same database.

I have tried arcpy.CopyFeatures_management and get the following error:

  ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000732: Input Features: Dataset \Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.dbo.AF_VW_CCTVwastewater does not exist or is not supported
  WARNING 000725: Output Feature Class: Dataset \Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.sde\TaupoGISTest.DBO.CCTV\TaupoGISTest.DBO.CCTV_Polyline_Copy already exists.
  Failed to execute (CopyFeatures).

I have tried arcpy.Append_management and get the following error:

  ExecuteError: Failed to execute. Parameters are not valid.
  ERROR 000732: Target Dataset: Dataset polyline does not exist or is not supported
  Failed to execute (Append).

I have tried arcpy.CopyRows_management and get the following error:

  ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Input Rows: Dataset
  \Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.dbo.AF_VW_CCTVwastewater
  does not exist or is not supported WARNING 000725: Output Table:
  Dataset
  \Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.sde\TaupoGISTest.DBO.CCTV\TaupoGISTest.DBO.CCTV_Polyline_Copy
  already exists. Failed to execute (CopyRows).

This is my code:

import arcpy
import os

arcpy.env.overwriteOutput = True

source_path = r'\\Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.dbo.'
output_path = r'\\Virgis01\gis\Naomi-AssetFindaTest\TaupoGISTest.sde\TaupoGISTest.DBO.CCTV\TaupoGISTest.DBO.'

polyline = output_path+r'CCTV_Polyline_Copy'

views = [r'AF_VW_CCTVwastewater', r'AF_VW_CCTVstormwater']

for view in views:
    print 'processing ' + view + '             o   o'
    print 'deleting                                      >  '
    #arcpy.TruncateTable_management(polyline)
    print 'copying                                    \_____/'
    #arcpy.CopyRows_management(source_path+view, polyline)
    #arcpy.Append_management("source_path+view", "polyline", "TEST")
    arcpy.CopyFeatures_management(source_path+view, polyline)

Are there another alternative that will turn this query layer into a feature class and bypass this error?


texts: Leaflet.contextmenu plugin runs into Uncaught Type error under certain conditions

texts:
I'm working on a leaflet map application getting data from OpenStreetMap and use the Leaflet plugin Leaflet.contextmenu.
I have defined a context menu on the map:

var Map = L.map('Map', {
    zoomControl: false,
    attributionControl: false,
    contextmenu: true,
    contextmenuWidth: 140,
    contextmenuItems: [
        {
            text: 'New ...',
            icon: 'img/pencil.png',
            retinaIcon: 'img/pencil_@2.png',
            callback: newFeatureOnMap
        }, {
            separator: true,
        }, {
            text: 'Address',
            icon: 'img/address.png',
            retinaIcon: 'img/address_@2.png',
            callback: showAddress
        }]
}).setView([51.0, 9.9], 6);

Additionally I declare a second context menu on markers on a L.geoJSON layer definition:

function createLayer (geoJsonData, networkName) {
    return L.geoJSON(geoJsonData, {

        filter: function (feature, layer) {
            return (checkAtmNetwork(networkName, feature));
        },

        pointToLayer: function (feature, latlng) {

            return L.marker(latlng, {
                    icon: L.AwesomeMarkers.icon(
                        {
                            icon: feature.properties.tags.amenity === 'bank' ? 'institution' : 'euro',
                            prefix: 'fa',
                            markerColor: MARKER_COLOR[networkName]
                        }),
                    contextmenu: true,
                    contextmenuInheritItems: false,
                    contextmenuItems: [{
                        text: 'Edit ...',
                        icon: 'img/pencil.png',
                        retinaIcon: 'img/pencil_@2.png',
                        callback: editFeature
                    }, {
                        text: 'Delete',
                        icon: 'img/trash.png',
                        retinaIcon: 'img/trash_@2.png',
                        callback: deleteFeature
                    },
                        {
                        separator: true
                    }, {
                        text: 'Route',
                        icon: 'img/route.png',
                        retinaIcon: 'img/route_@2.png'
                        callback: routeToFeature
                    }]
                }
            );

        },
   });

Both context menus work independent as they should do. But, if I move the map quite often (on moveend event the layers a completely removed from the map and newly loaded and created), sometimes the context menu of the markers are doubled in the DOM and a right click on marker produces the error:

leaflet.contextmenu.js:548 Uncaught TypeError: Cannot read property 'contextmenu' of null
    at e._hideContextMenu (leaflet.contextmenu.js:548)
    at e.fire (leaflet.js:5)
    at e._hide (leaflet.contextmenu.js:393)
    at HTMLAnchorElement.<anonymous> (leaflet.contextmenu.js:324)
    at HTMLAnchorElement.h (leaflet.js:6)

The context menu on markers won't work any more, and on the map is shown the marker context menu, but not working too. The map context menu is not shown any more.

May be, there are problems with the life-cycle of the marker layer on moving the map and reloading the data...

Has anyone an idea what's going wrong?


texts: Speed optimization of shadow computation with grass gis r.sun

texts:
I am trying to optimize the speed of r.sun to compute shadowed area on a DTM at 1m resolution. My current area is 3.39 km² but my final area will be 450 km².

The r.sun wiki say that pre-computing horizons as well as lat/lon increase the model, but the opposite happen on my computer (GRASS 7.0.1RC2, Windows7 64 bits [although r.sun run in 32 bits mode], 4 CPU, 8 threads)

Here is an example:

r.sun elevation=MNEcomplet@ShadowProject day=88 time=7 incidout=incidoutD88T07 --overwrite npartitions=2 distance_step=1

==> 38 sec

r.sun elevation=MNEcomplet@ShadowProject horizon_basename=horizon horizon_step=30 day=88 time=7 incidout=incidoutD88T07 --overwrite npartitions=2 distance_step=1.0 lat=MNEcomplet.lat long=MNEcomplet.lon

==> between 56 sec and 1min7s

npartitions are different on the 2 examples (I have no choice), but this parameter has not a big impact on the computing time.

Does anyone have an idea why it has this behavior?

Testing on Grass 6.4, using horizons allows me to decrease computation time from 42-50s to 32s.

Another question is: Is there a possibility to avoid the computation of the sun angles for each cell ?

Providing a single value for the whole area would be ok for me, and I can compute it by myself


texts: Filling contour levels with grain (texture) using Leaflet

texts:
I created a program to draw temperature contour map of France:

This is the code:

library(ggplot2)
library(gstat)
library(sp)
library(maptools)
require(rgeos)
library(rgdal)
library(foreign)
library(maptools)
library(knitr)
library(jsonlite)
library(raster)
library(leaflet)

#récupérer les données a partir de json
dataJson <- fromJSON(readLines('jsonallcoor.js'))
# Téléchargement du fond de carte
FRA=readShapePoly("FRA_adm0.shp")

#concour (interpolation idw)
x=data.frame(0)
y=data.frame(0)
t=data.frame(0)
for ( i in 1:103) {
  x[i]<- as.data.frame.factor(dataJson[[i]]$longitude)# define x & y as longitude and latitude
  y[i]<- as.data.frame.factor(dataJson[[i]]$latitude)
  t[i]<- as.data.frame.factor(dataJson[[i]]$temperature)
}
#changer le nom de dataframe

x<- as.numeric(x)
y <- as.numeric(y)
temperature<- as.numeric(t)
frame=as.data.frame(cbind(x,y,temperature))
frame.xy = frame[c("x", "y")]
coordinates(frame.xy ) <- ~x+y
plot(frame.xy )

#Define the grid extent
x.range <- as.numeric(c(-4.445833,9.484722))  # min/max longitude of the interpolation area
y.range <- as.numeric(c(40.50306,51.2))  # min/max latitude of the interpolation area

#Create a data frame from all combinations of the supplied vectors or factors. See the description of the return value for precise details of the way this is done. Set spatial coordinates to create a Spatial object. Assign gridded structure:
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 0.1),
                   y = seq(from = y.range[1], to = y.range[2], by = 0.1))
# expand points to grid
coordinates(grd) <- ~x + y
#to pixel
gridded(grd) <- TRUE
#Plot the weather station locations and interpolation grid:
plot(grd, cex = 1.5, col = "grey")

#les lon et lat
points(frame.xy, pch =1, col = "black", cex = 0.1)
#Interpolate surface and fix the output
idw <- idw(formula = temperature ~ 1, locations = frame.xy,
           newdata = grd)  # apply idw model for the data

residual_grid = raster(idw, "var1.pred")
contourr <- rasterToContour(residual_grid)

library(leaflet)
## Initialisation
m <- leaflet(padding = 0)
m <- addTiles(m)## Ajout des pays-
m <- addPolygons(map = m, data = FRA, opacity = 100,
                 color = "#FAFCFA",
                 weight = 0.25,popup = NULL,
                 options = list(clickable = FALSE),
                 fill = T, fillColor = "#B3C4B3",
                 fillOpacity = 100)
## Ajout des cercles
for ( i in 1:103) {
  m <-addCircleMarkers(map = m,
                       lng = dataJson[[i]]$longitude,
                       lat = dataJson[[i]]$latitude,
                       radius=5, weight = 0.25,
                       stroke = T, opacity = 100,
                       fill = T, fillColor = "#920000",
                       fillOpacity = 100,
                       popup = dataJson[[i]]$station,
                       color = "white")
}
m <- fitBounds(map = m,
               lng1 = -4.445833,
               lat1 = 40.50306,
               lng2 = 9.484722,
               lat2 =  50.56417)
m<-addPolylines(map=m,data =contourr,fillOpacity=5,fillColor = "transparent",opacity=10,weight=1)
m<-addRectangles(map=m,
                 lng1=9.484722, lat1=40.50306,
                 lng2=-4.445833, lat2=51.2,
                 fillColor = "transparent"
)

## Dimensions de la carte
m$width <- 1000
m$height <-700

library(htmlwidgets)
saveWidget(m, 'contour.html', selfcontained = FALSE)

This is the result:

Now, I am looking a method for filling the contour levels with grain (texture).

EDIT:

I show you what I've done
I added a grid raster then I intersected with the card

but my ultimate goal is to combine texture with temperature value so the change in the value changes with grain like this :

I don't know how do it. please give me suggestion


texts: Trying to identify where is THIS St. John's Way

texts:
I'm trying to identify where is THIS St. John's Way.

Can't find it with Google Maps, too many irrelevant places where showing up.


texts: Distance to nearest point for every point same SpatialPointsDataFrame in R

texts: I have a SpatialPointsDataFrame that I'm working with in R. I want to add a new vector to the dataframe that contains, for every point, the distance to the nearest other point in the SpatialPointsDataFrame. I've looked at knearneigh in the spdep package as well as spDistsN1 and spDists in the sp package, but neither of those give me exactly what I want. To be clear, all of these points are in ONE SpatialPointsDataFrame.


texts: Computing correlation between two image collections in GEE

texts:
I have two image collections with 744 images each (and the same resolution and extent), covering the whole Earth.
I need to compute the correlation between the image collections for each pixel.


texts:
I'm using worldview imagery downloaded from gbdx and I want to apply cloud mask for it. It's preferable to use ready-made layer or some python module/tool to create it. The demand for the method precision: it's allowable to cut clouded areas with some reserve (lose a part of useful information).

What should I try to consider?

I thought about some threshold, but perhaps it would be a bad idea.


texts: Automating polygon "useable area" check in QGIS?

texts:
I need to run a query that will identify all property remainders (available land for construction purposes) that can accommodate a 40m² building footprint (must be square).

I have used imagery classification to create building footprints for existing structures on properties and have clipped them from the land parcels for calculating area of undeveloped land within each property boundary.

I can't assume that each remainder measuring >40m² can accommodate the construction of a 40m² building because the calculation will include any sliver of land that wasn't part of the existing building area, for example, the gap between house and boundary, small front garden area etc.

I need to do this check for approximately 20,000 property boundaries, so automation is key and I have little to no coding experience.


texts: How to change the view of a map window in MapInfo to a new center coordinate?

texts:
How can I change the map view of a map window in MapInfo Professional to a new center coordinate? I need to change the viewed map in the MapInfo window via a predefined set of coordinates via MapBasic.

I tried to find it out through the MapBasic window but didn't find anything useful.


texts: Configure GDAL to use Open JPEG 2000 library

texts:
gdalinfo --formats shows one entry for the JPEG 2000:

JPEG2000 -raster,vector- (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1), based on Jasper library

I don't want to use the Jasper library. I want to use the Open JPEG 2000 library. How do I configure GDAL to use a different library for a particular raster format?


texts:
I run a command line gdal_calc.py in a scala script. I run this script with 'sbt run' on my terminal.
I got this error :

Traceback (most recent call last):
  File "/usr/local/bin/gdal_calc.py", line 329, in <module>
    main()
  File "/usr/local/bin/gdal_calc.py", line 326, in main
    doit(opts, args)
  File "/usr/local/bin/gdal_calc.py", line 282, in doit
    myResult = ((1*(myNDVs==0))*myResult) + (myOutNDV*myNDVs)
TypeError: only integer arrays with one element can be converted to an index

Running the exact same command line directly in my terminal works well.
The environment is the same : the gdal library used is the same

The command line runned is the following :

gdal_calc.py --outfile=outfile.tiff -A infile.tiff --overwrite --calc="3*(A==2)"

Can someone explain it ?


texts:
We are using a intranet geoloket (html). I would like to share some QGIS projects using PHP.
This is the code I start with:
header("Content-disposition: attachment; filename=MapInQgis.qgz");

header("Content-type: application/**octet-stream**");

readfile("MapInQgis.qgz");

Which content-type must I use to make sure my collegues can open the MapInQgis.qgz-file?


texts: Performing Geographic Weighted Regression (GWR) in Java?

texts:
Is there any existing GWR (Geographic Weighted Regression) packages in Java?

The packages should be able to conduct the basic GWR, the bandwidth optimation and some significance test.


texts: Get the name of the layer from a layer Vector

texts:
My aim is to highlight that two differently colored layers wfs when passing the mouse over one of these layers.

It would need to obtain the name of the Vector Layer and also use a style in one layer and another in another.

This is my code WFS:

  var Urbana_source = new ol.source.Vector({
             format: new ol.format.GeoJSON(),
             url: function(extent, resolution, projection) {
      return 'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=pucol:Urbana&' +
    'outputFormat=application/json&srsname=EPSG:3857&' +
    'bbox=' + extent.join(',') + ',EPSG:3857';
  },
  strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
  maxZoom: 19
  }))

   });
   var urbana_wfs = new ol.layer.Vector({
                    source: Urbana_source,
                    style: new ol.style.Style({

                   stroke: new ol.style.Stroke({
                  color: 'rgba(0, 0, 0, 0)',
                  width: 2
       })
     })
   });
   urbana_wfs.set('name', 'urbana');

  var Rustica_source = new ol.source.Vector({
  format: new ol.format.GeoJSON(),
  url: function(extent, resolution, projection) {
   return 'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=pucol:Rustica&' +
    'outputFormat=application/json&srsname=EPSG:3857&' +
    'bbox=' + extent.join(',') + ',EPSG:3857';
   },
    strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
    maxZoom: 19
   }))

 });
 var rustica_wfs = new ol.layer.Vector({
 source: Rustica_source,
 style: new ol.style.Style({

  stroke: new ol.style.Stroke({
     color: 'rgba(0, 0, 0, 0)',
        width: 2
     })
   })
 });
   rustica_wfs.set('name', 'rustica');

    var displayFeatureInfo = function(pixel) {

   var feature = map.forEachFeatureAtPixel(pixel, function(feature, layer) {
    return feature;
  });

   var info = document.getElementById('info');
   if (feature) {
  //info.innerHTML =feature.get('referencia');
    var ref =feature.get('referencia');
   document.getElementById("Getfeatureinfo").value = ref;
   document.getElementsByName("info").innerHTML=ref;
  } else {
  info.innerHTML = '&nbsp;';
 }

 if (feature !== highlight) {
   if (highlight) {
    featureOverlay.getSource().removeFeature(highlight);
  }
  if (feature) {

    featureOverlay.getSource().addFeature(feature);
  }
map.getLayers().forEach(function(el) {
      if (el.get('name') === 'urbana') {
        console.log(urbana_wfs.get('name'));
      }else if (el.get('name') === 'rustica') {
        console.log(rustica_wfs.get('name'));
      }
    })
highlight = feature;
  }

  };
  map.on('pointermove', function(evt) {
 if (evt.dragging) {
 return;
 }
 var pixel = map.getEventPixel(evt.originalEvent);
 displayFeatureInfo(pixel);
  });

 map.on('click', function(evt) {
 displayFeatureInfo(evt.pixel);
 });
 };

This is my code.
What I would like to make is that when I pass the mouse over the layer urbana_wfs color glows and when I do the other rustic lighting.

My problem is that when I do this, I always get the two log, not one.
as you could get the name of this layer?


texts:
I saw a great interface to FME with Python

What you guys are doing with it? I want ideas.


texts: Append daily rasters to monthly lists in Python Arcpy

texts:
I'm trying to go over a folder of geotiff rasters where the last three digits of the raster name is the julian date and I want to append the correct days to monthly lists (Jan-Dec) and then run Cell Statistics on each monthly list of rasters. But I can't get the rasters to append to the monthly lists in Python.
I'm using ArcGIS 10.3.1 and Python 2.7.8.

What I have so far isn't working:

import arcpy
from arcpy import sa
from arcpy import env
import glob

directory = myFilePath

lstrasters = glob.glob(directory + os.sep + "*.tif")
lstrasters.sort()
print lstrasters

m1 = [] #empty Jan list
m2 = [] #empty Feb list
m3 = []
m4 = []
m5 = []
m6 = []
m7 = []
m8 = []
m9 = []
m10 = []
m11 = []
m12 = []

for lstras in lstrasters:
    jdate = str(lstras[-7:-4])
    day = int(jdate)
    print day
    sa.Raster(lstras)
    if day in range(1,32):
        m1.append(lstras)
    if day in range(32,60):
        m2.append(lstras)
    if day in range(60,91):
        m3.append(lstras)
    if day in range(91,121):
        m4.append(lstras)
    if day in range(121,152):
        m5.append(lstras)
    if day in range(152,182):
        m6.append(lstras)
    if day in range(182,213):
        m7.append(lstras)
    if day in range(213,244):
        m8.append(lstras)
    if day in range(244,274):
        m9.append(lstras)
    if day in range(274,305):
        m10.append(lstras)
    if day in range(305,335):
        m11.append(lstras)
    if day in range(335,366):
        m12.append(lstras)
    else:
        print "jdate is out of range"

print m1
print m2

....etc.

It prints out the correct jdate day for each raster for sure. But I can't get any of the rasters to append to any of the empty monthly lists that match up to the jdate. it just prints the 'else' statement for each raster saying it is out of range and then prints out each m1,m2,... as empty lists.

What am I doing wrong?

Edit: I changed the day = jdate to day = int(jdate)

Edit 2: Based on the comments, I made the following change to the if/else statements and the script is working properly now.

for lstras in lstrasters:
    jdate = str(lstras[-7:-4])
    day = int(jdate)
    print day
    sa.Raster(lstras)
    if day in range(1,32):
        m1.append(lstras)
    elif day in range(32,60):
        m2.append(lstras)
    elif day in range(60,91):
        m3.append(lstras)
    elif day in range(91,121):
        m4.append(lstras)
    elif day in range(121,152):
        m5.append(lstras)
    elif day in range(152,182):
        m6.append(lstras)
    elif day in range(182,213):
        m7.append(lstras)
    elif day in range(213,244):
        m8.append(lstras)
    elif day in range(244,274):
        m9.append(lstras)
    elif day in range(274,305):
        m10.append(lstras)
    elif day in range(305,335):
        m11.append(lstras)
    elif day in range(335,366):
        m12.append(lstras)
    else:
        print "jdate is out of range"


texts: CartoDb Response for preflight is invalid (redirect) error

texts:
I created a map in carto and then applied to my c# application like CartoDB source example that I saw on openlayers web page.  like this:

     var mapConfig = {
        'layers': [{
            'type': 'cartodb',
            'options': {
                'cartocss_version': '2.1.1',
                'cartocss': '#layer { polygon-fill: #374C70; polygon-opacity: 0.9; ::outline { line-color: #FFF; line width: 1; line-opacity: 0.5; } }',
                'sql': 'SELECT * FROM kitap1'
            }
        }]
    };

    var cartoDBSource = new ol.source.CartoDB({
        account: 'seymacengiz',
        config: mapConfig
    });

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

At first, it was working .But now I'm getting this error:

    XMLHttpRequest cannot load https://seymacengiz.cartodb.com/api/v1/map. Response for preflight is invalid (redirect)

And also this link goes to "https://carto.com/noneuser.html"
I saw v1 in link.I think, it means version1 .

Can this be the reason for the error?


texts: POSTGIS raster/Java with GDAL bindings - Load byte array in Gdal Dataset Object

texts:
I would like to load a raster as  java.awt.image.BufferedImage (or equivalent class that would allow to perform some display) out of my PostGIS database. My Raster is inDb (loaded with raster2pgsql)

I already created a proper DAO to "connect" the Postgis database and my Java package.
So far, I managed to import an display tiles from my raster using POSTGIS ST_AsPNG queries as presented in the code extract below.

public BufferedImage test_raster() {

    byte[] content = null;
    Dataset dattaset = null;

    try {

        ResultSet result = this .connect
                                .createStatement(
                                            ResultSet.TYPE_SCROLL_INSENSITIVE,
                                                ResultSet.CONCUR_UPDATABLE
                                         ).executeQuery(
             "SELECT ST_AsPNG(rast) as tmp FROM bath50m  WHERE ST_Intersects(rast, ST_GeomFromText('MULTIPOLYGON(((50 50, 50 0, 0 0, 0 50, 50 50)))',4326))"
                                         );
        while(result.next()){
            if (result.absolute(result.getRow())){

            content = result.getBytes("tmp");
                ByteArrayInputStream bis = new ByteArrayInputStream(content);
            BufferedImage image = ImageIO.read(bis);

            }
        }

    } catch (SQLException e) {
                e.printStackTrace();
    }
    return image;
}

That is interresting but I'd like to keep the geographic metadata that are embedded with the raster.

I decided to give a try to extract tiles using the POSTGIS ST_AsGdalRaster so metadata would be passed along with the image. I properly compile GDAL to get the java bindings and set PATH so that I can use it (I tried succesfully the Nasa Worldwind GDALtest.java)

According to the doc, it returns the raster tile in the designated GDAL Raster format. Raster formats are one of those supported by your compiled library. Use ST_GDALRasters() to get a list of formats supported by your library.

Actually ST_AsGdalRaster returns a bytes array that I can't cast to the Gdal Dataset type. My intention was to adapt the openfile method from the GDALtest.java  to pass byte array or whatever I could.

I am not sure where I am going with that. I know I am doing something wrong but I miss support/experience to realize it.

Any help would be welcome.

Thanks for reading me

edit: I realize I need to clarify meaning. I try to use the GDALDataset Class with a byte array as input. Is it possible?


texts: Editing joined parcels to eliminate gaps in parcel fabric?

texts:
I have a parcel fabric layer that has a gap in it:

I want to edit the selected parcel on the bottom to fill the gap. I don't want to construct a new parcel to do so. How do I modify the parcel in such a way that it maintains its attributes and fills in the gap with the parcel above it?

Based on answers and comments, I created a topology and placed the parcel fabric inside it. It appears that you cannot use topological tools to edit a parcel fabric. As such, any answers using topology will not work.


texts: QGIS Field Calculator Aggregation with a intersect query

texts:
I have point (manholes) and line (sewer main) datasets and want attributes from the points at both ends of the line.

If I run the following after I install RefFunctions in field calc on the line

geomintersects('Sewer_MH_Test','RL_of_Lid')

I get the value I need from the start of the line (3.78). How do I get the value at the end (3.75) into a new column? Is there some aggregator to use?

Based on @Taras' answer I have been trying the following but it keeps crashing.

WITH start_data AS (
SELECT l.GID_New AS line_id_start, p.US_IL AS x
FROM points AS p
JOIN lines AS l ON st_intersects(start_point(l.geometry),p.geometry)
),

end_data AS (
SELECT l.GID_New AS line_id_end, p.DS_IL AS y
FROM points AS p
JOIN lines AS l ON st_intersects(end_point(l.geometry),p.geometry)
)

SELECT l.*, start_data.x, end_data.y
FROM lines AS l
JOIN start_data ON l.GID_New = line_id_start
JOIN end_data ON l.GID_New=line_id_end

This shows the output I want...

The test data is only about 20 lines and points - it can be downloaded from https://drive.google.com/file/d/12rUV_pAaevs9vATbIsgc9QJvvziRygGN/view?usp=sharing


texts: Determine acquisition date from a QGIS ESRI Basemap

texts: I am using the ESRI Satellite base map on QGIS. I need to know its acquisition date, is it possible? And if yes, how?


texts: Is there a way to use a different symbol on a vertical line vs a horizontal line in QGIS?

texts:
I'm working on a map that has an oldtime feel. There are cliffs on the map, some of them are roughly horizontal (running east to west) and some of them roughly vertical (running north to south) that I have indicated using a line. I would like to use a different symbol on the vertical line segments than I use on the horizontal line segments.
I would create another layer, one with horizontal and one with vertical line segments, but this is inconvenient because many of the cliffs curve around. I was hoping there would be a more simple method.
Thanks for the help!
Edit: here's a quick sketch of what I'm trying to accomplish via lines


texts:
I'm trying to add in QGIS a layer shown here : propluvia

When I watch in my console I see this request corresponding to the layer I need.

The server is not able to respond to a getCapabilities request for whatever reason.

So I tried to load this layer via this GDAL XML file.

<GDAL_WMS>
<Service name="WMS">
    <ServerUrl>http://eau.agriculture.gouv.fr/WMS-T.php?map=/export/home/rimboval/geomap/mapfiles/propluvia.map&</ServerUrl>
    <SRS>EPSG:2154</SRS>
    <ImageFormat>image/png</ImageFormat>
    <Transparent>TRUE</Transparent>
    <Layers>Zones_Niveau_Sup</Layers>
    <Change key="${date}">2014-08-18</Change>
</Service>
<DataWindow>
    <UpperLeftX>508204</UpperLeftX>
    <UpperLeftY>6363402</UpperLeftY>
    <LowerRightX>627904</LowerRightX>
    <LowerRightY>6468402</LowerRightY>
    <SizeX>256</SizeX>
    <SizeY>256</SizeY>
</DataWindow>

I'm using key="${date}" because I saw it here.

The resulting image looks like in this question.

I tried to watch what QGIS is requesting with fiddler but I can't see any request while it is working with a valid WMS layer.


texts: Using attribute domain and subtypes in ArcMap

texts:
I have a schema (attached) for a certain classification that I need to include in the data structure.
Is it possible to use subtypes and domains to get the proper solution?
Simply, How to represent schema levels in database structure with a minimum number of fields?

I am using ArcMap 10.4.


texts: Scale raster in ArcGIS using raster calculator

texts:
I have a raster with values ranging from 0 to 8991130624. This raster is actually the result of a cost distance analysis.

I want to scale the raster to range from 0 to 1 and take the inverse? so that a cell with the highest cost/value (e.g. 8991130624) will now have a value of 0 while a cell with the lowest cost/value (e.g. 0) will now have a value of 1.

I am using raster calculator with the following expression:

1-(x/8991130624)

where I add x to the expression by clicking on the raster from the list of available rasters when using the tool

This is not working and producing an error:
"ERROR 000539: Error running expression; rcexec() :ERROR 999998:Unexpected Error."

Can someone tell me what I'm doing wrong?


texts: I wanted to install latest version of GDAL (2.3) in Python 2.7.15 32 bits. I usually install the gdal core and bindings from www.gisinternals.com for my python version (msc 1500). However, for version of GDAL 2.3 there is only available for msc 1900 and 1911 which I dont know what version of Python correspond to, but not mine (maybe Python 3?). So, is there no GDAL 2.3 for version 2.7.15? Anyone knows?


texts: How do I convert to EPSG:3857 with GDAL?

texts:
I have some TIFs in EPSG4326 the GeoTiff info looks like:

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]

How do convert them into EPSG:3857 so that the GeoTiff info reads something like:

PROJCS["Google Maps Global Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]


texts: Can conversion of MapInfo WOR to ArcGIS MXD preserve styles and coordinate systems?

texts:
Is there any ways to migrate from MapInfo's data with WOR project file to ArcGIS MXD project file without leaving all styles and coordinate system?
I know that I can convert from TAB files to SHP files via ogr2ogr but what's about the styles information?

P.S. I have MapInfo Professional 12 and ArcGIS Desktop 9.3.1. Also, I can't use the Tab Reader Tool now.


texts: Cropping GeoTIFF with GeoTools - remove NO_DATA border

texts:
I have a raster with borders full of NO_DATA value.
Is it possible to crop a raster over a "minimum" bounding rectangle,

Legend: * - no data, + - data

*******
* ++**   ->  ++
* ++**       ++
*******


texts:
I wrote this script to classify a shapefile in 5 classes, this works for small files (minor 8000 features), for larger files, the process has stopped. Can make this better?

#Load the shape
vrAdubo = QgsVectorLayer(vetor_saida_adubo, "Adubo", "ogr")

if not vrAdubo:
   progress.setText("Erro ao carregar shape resultante Adubo")

#Insert fields area_ha,  classe e kmph
provVrAdub = vrAdubo.dataProvider()
provVrAdub.addAttributes([QgsField("area_ha", QVariant.Double), QgsField("classe", QVariant.Int), QgsField("kmph", QVariant.Double), QgsField("area_class", QVariant.Double),
QgsField("perc_class", QVariant.Double)])
vrAdubo.updateFields()

#Set fields valures
vrAdubo.startEditing()
iterVrAdub = vrAdubo.getFeatures()

#variaveis areas totais por classe
area_cls1 = 0.000
area_cls2 = 0.000
area_cls3 = 0.000
area_cls4 = 0.000
area_cls5 = 0.000

#Setting filds
for featAdubo in iterVrAdub:
    geom = featAdubo.geometry()
    featAdubo['area_ha'] = geom.area() / 10000
    areatt_ad += featAdubo['area_ha'] #area total mapeada
    ct_pol_adu += 1 #contagem de poligonos
    featAdubo['kmph'] = featAdubo['Speed'] * 3.6

#Classify by ranges
    if featAdubo['AppldRate'] > 0.0000:
        if featAdubo['AppldRate'] < 400.000:
            featAdubo['classe'] = 1
            area_cls1 += featAdubo['area_ha']
        elif featAdubo['AppldRate'] >= 400.000 and featAdubo['AppldRate']< 450.000:
            featAdubo['classe'] = 2
            area_cls2 += featAdubo['area_ha']
        elif featAdubo['AppldRate'] >= 450.000  and featAdubo['AppldRate']< 550.000:
            featAdubo['classe'] = 3
            area_cls3 += featAdubo['area_ha']
        elif featAdubo['AppldRate'] >= 550.000  and featAdubo['AppldRate']< 600.000:
            featAdubo['classe'] = 4
            area_cls4 += featAdubo['area_ha']
        else:
            featAdubo['classe'] = 5
            area_cls5 += featAdubo['area_ha']
    else:
        progress.setText("Aplicacao zerada em AppIdRate: %" % featAdubo.id())

    vrAdubo.updateFeature(featAdubo)

vrAdubo.commitChanges()


texts: How to load a style from postgres on a sub request in PyQt

texts:
In a qgis-processing script, I'm loading sub-requested layers with x, y, radius to have less data in qgis-desktop. The data is well loaded but I can't apply styles stored in postgresql layer_styles table.
The Qgis layers'name, the f_table_name and the stylename are the sames. Could somebody help me finishing this script?

##myscript=group
##x=number 449000
##y=number 219000
##radius=number 500

from qgis.core import *
from PyQt4.QtCore import *

# list of layers' names
file = open(r"\\fileserver\sig$\symbols\mapping.csv", "r")

for line in file:
    registry = QgsMapLayerRegistry.instance().mapLayers()

    # remove existing layers
    for layer in registry.values():
        if (layer.name() == line):
            QgsMapLayerRegistry.instance().removeMapLayer(layer)

    # making the sub request
    part1="(select * from schema."
    part2=" WHERE ST_INTERSECTS(ST_SetSRID(ST_Buffer(ST_MakePoint("
    part3=","
    part4="),"
    part5="),3163),geom))"
    sql = part1 + line + part2 +str(x)+ part3 +str(y)+part4+str(radius)+part5

    # building the uri
    uri = QgsDataSourceURI()
    uri.setConnection("192.168.1.71","5432","mybase","me","mypass")
    uri.setDataSource("",sql,"geom","","gid")

    # adding the layer to the registry
    vlayer = QgsVectorLayer(uri.uri(),line,"postgres")
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)
    # applying the style
    vlayer.applyNamedStyle(line)

    if not vlayer.isValid():
        print "layer failed to load"


texts: How to covert Gauss-Kruger coordinates to UTM coordinates with GeoTools java library? As I have figured out the location (zone, country) is important. I need the calculation for Slovenia.


texts: How do I add vertices on a segment in QGIS 3.9.? Double-clicking doesn't work

texts: I tried adding vertices to segments of a polygon by double-clicking on the segment, but instead of adding a node, the nearest node was moved to the double-click. Clicking on the cross doesn't work either as this is only visible in longer segments.


texts: How do I insert new-line character into Addin Extension Description?

texts:
Apparently the text box in Esri's Extension manager doesn't wrap text (10.0 SP2).

As a workaround, I'm trying to insert newline characters into the xml file, but haven't found anything that works.

Here's the xml:

<Extensions>
  <Extension id="AmberGIS_ExtensionMgrAddin_ExtMgrExtension" class="ExtMgrExtension" productName="Extension State Manager" showInExtensionDialog="true" autoLoad="true">
    <Description>Lets you specify extensions that will be \n turned off automatically when document closes.</Description>
  </Extension>
</Extensions>

Does anyone know how to put a newline character into the add-in xml config file?


texts: I want to Split by Attribute to run a function multiple times, but 'Id' is all 0?

texts:
I want to run my function over each feature in a point shapefile layer, but splitting by attribute returns the same number of points, because the Id field is all 0.
Is there a way, using ArcPy syntax, to add unique values to the Id field, or an entirely new field of unique values?
Or, a way to split by the FID attribute, which is always unique values?
    import arcpy
    from arcpy import env
    from arcpy.sa import *
    #DEFINE PARAMETERS & zLimit
    env.workspace = arcpy.GetParameterAsText(0)     #type is folder
    inSurfaceRaster = arcpy.GetParameterAsText(1)   #type is Raster Layer
    CROSSINGS = arcpy.GetParameterAsText(2)         #type is Table View
    EXTENT = arcpy.GetParameterAsText(3)            #type is double
    PARCEL = arcpy.GetParameterAsText(4)            #type is Table View
    zLimit = 3.2808
    #CHECK OUT REQUIRED EXTENSION
    arcpy.CheckOutExtension("Spatial")
    #CLIPPING RASTER LAYER TO SMALLER AREA
    AOI = arcpy.Buffer_analysis(PARCEL,'AOI',10000)
    USR = ExtractByMask(inSurfaceRaster,AOI)
    #GENERATING REQUIRED RASTER LAYERS
    outFill = Fill(USR,zLimit)
    outFlowDir = FlowDirection(outFill)
    outFlowAcc = FlowAccumulation(outFlowDir)
    outSnapPoint = SnapPourPoint(CROSSINGS, outFlowAcc, EXTENT)
    #DELINEATION OF WATERSHEDS
    watersheds = Watershed(outFlowDir,outSnapPoint,"VALUE")
    arcpy.RasterToPolygon_conversion(watersheds,"BASINS")
    arcpy.AddMessage("BASINS"+CROSSINGS)
    #DELETE EXCESS SHAPEFILES
    arcpy.Delete_management(outFill)
    arcpy.Delete_management(outFlowDir)
    arcpy.Delete_management(outFlowAcc)
    arcpy.Delete_management(outSnapPoint)
    arcpy.Delete_management(watersheds)
    arcpy.Delete_management(AOI)
    arcpy.Delete_management(USR)

The above code draws a drainage basin for the 'CROSSINGS' input, my issue is that only one basin is drawn, when the idea is to have one for each point feature in 'CROSSINGS'.
If I could split 'CROSSINGS' I would then run the function over each new point.


texts: Draw line drawn from point to an open infowindow?

texts:
I want the infowindows on my map to be in a fixed location, but that sometimes makes it hard to tell which point was clicked (which point is associated with the infowindow displayed).
Is there a way to draw a line from the point to an opened infowindow located in a fixed position or make the selected point stand out from the rest of the non-selected points?


texts: Plotting Geographical Grid from Long/Lat Coordinates in R

texts:
I have a file of coordinates (Long/Lat) that should enable me to plot a series of 10km grid squares covering the UK. However, for each grid square I have five sets of coordinates by way of plotting the four corners of the square (starting bottom left and moving round clockwise) and returning to the starting point (so the first and last sets of coordinates for each square are identical), as follows:

> head(Coordinates)
       long      lat order grid resolution

1 -2.001839 60.68321     1 HP00         10
2 -2.001845 60.77299     2 HP00         10
3 -1.818269 60.77287     3 HP00         10
4 -1.818774 60.68309     4 HP00         10
5 -2.001839 60.68321     5 HP00         10
6 -2.001845 60.77299     1 HP01         10

I've been exploring options for plotting the grid (e.g. geom_rect, geom_raster etc.), but the fact that I have five sets of coordinates for each square ('order' in the table above) seems to be complicating matters! That said, I'm very new to R, so may well be missing something that's obvious to others. Can anyone suggest how I might best go about using the data to generate/plot the complete set of grid squares, please?


texts: How to create a valid global polygon grid in PostGIS?

texts:
I want to create a grid spawning the whole globe. The grid is used to measure a density value by counting the number of certain spatial objects in each grid cell.

To create the grid, I have created multiline polygons, both using a script of my own, and the function suggested here: How to create a regular polygon grid in PostGIS? . Both algorithms loop from -90 .. 90 , -180 .. 180 in a fixed step size creating the multiline polygons.

I have two problems.

First, adjacent cells share a separating line, making objects on it counting towards two or more cells. I could subtract a very small value from corners, but this risks loosing objects when they come to lay exactly on the original line. The question is thus if there is a constant for the absolute minimum step value that I could use to guarantee that one point exactly lays next to another with no space inbetween.

Second, it seems the algorithm generates invalid geometries. e.g when I compute the square meters of a grid (for the density)

SELECT ST_Area(ST_GeographyFromText('POLYGON((-15 0,-15 5,-10 5,-10 0,-15 0))'))

throws an error ptarray_area_spheroid: cannot handle ptarray that crosses equator
If I avoid the geography type, and use a geometry there is no global projection that outputs the correct size of all grids. e.g

Select ST_Transform( ST_GeomFromEWKT('SRID=4326;POINT(-180 90 0)'), 32630 )

throws couldn't project point (-180 90 0): latitude or longitude exceeded limits (-14)

How could I generate a grid of valid geometries that covers the whole globe?


texts: I have been trying to warp ECW and Ozfx3 map formats to GeoTiff format.  I have navigated to Raster>Projections>Warp(reproject) I have then selected the file to warp the output that I want and the process starts but each time the process fails within 5  seconds with the dialogue  "The process crashed some time after starting successfully".


texts: Is there any comprehensive book or online instruction on image segmentation techniques?

texts:
I'm looking for a source to learn theoretical basis of image segmentation techniques thoroughly.

The desirable book, or other kinds of source you may know, should cover the techniques that are more applicable in remote sensing studies.


texts:
I want to give a map that contains hyperlinks as a package of Adobe or
an exe file to someone who does not have access to ArcMap. This
package should have just some points and polygons that when he click
on each of them, a file which is linked to that point or polygon
opens.

How can I do that?


texts:
I want to convert large map file to tms format.
how directory structure, zoom levels , tile size, xml file creation have to be dealt with.
Please point to the relevant docs.

And is there any gdal c++ api to create tiles(not gdal2tiles.py).
Thanks


texts:
Possible Duplicate:
  Are long lines in shapefiles to be considered geodesics or straight lines in the 2D latlong space?

In various representations, such as shapefiles and kml, polygons are defined by their vertices. My question is what is the nature of the "line" which links two vertices? Is it a great circle arc? Is it whatever is a straight line in the current projection? If the latter, then the actual interior of the polygon in a function of the current projection, which seems undesirable.


texts: Ordnance Survey Master Map font in QGIS 2.14

texts:
I have been trying to style some OSMM data using their supplied QML files which include a font with some of their symbology in. I've added the font to my system (ubuntu 14.04 LTS) and when viewed in FontMatrix and other system font tools it looks fine:

But when I come to use it in QGIS the font has no glyphs visible and I end up with ugly rectangles instead of my symbols.

Is there some trick to get QGIS to recognise this font?


texts:
I have two vector layers as seen below, one is a landcover map showing the different habitats in a region and one is a grid of 4 x 4 km squares that covers the region.

I have to extend this grid to cover the whole landcover layer as currently it does not cover it properly, I could also create a new one but I need the current squares to stay in the same place. I can't find anything on expanding this layer or to create a new one based on the old one


texts:
I have download raw landsat 8  image with Rasterio. The image was download from GEE. . The image has 11 bands.
I'm trying to display the image  with grey scale/color scale as DN values, but for some reason I get the image in white values.
This is how I open the image and try to dusplay it:
#import required libraries
import rasterio
from rasterio import plot
import matplotlib.pyplot as plt
from rasterio.plot import show
%matplotlib inline

#open the raster
src=rasterio.open('16092020_landsat.tif')

#display one band:
plt.imshow(src.read(1), cmap='pink')

#RGB image
plt.figure(figsize=(20,10))
show(src.read([4,3,2]),transform=src.transform,title='Image- bands 4,3,2 IN 0-255 Values')

I get this image with correct coordinates but is white.
When I try to plot the pixels as histogram with the RGB I get empty plot:

When I run the same script with sentinel2 reflectance all works so I don't understand why here I get nothing,
My end goal: to be able to show the Landsat8 image as well.
edit: I suspect it has to do with the values of the Landat as they are between 0-25k


texts: What meter resolution are US Census Block Shapes?

texts:
I'm creating a map layer that supports up to zoom level 14 and I want to know what resolution US Census blocks are?

At zoom level 14 each pixel is approximately 16 meters.

I am wondering if I can simplify the block data to 20m resolution or if it is worthwhile to go with a higher resolution?


texts: Reading column line by line from table in geodatabase using ArcPy?

texts:
I have a table in a geodatabase.

It looks like this:

   id   table_na
0   1   ER_ASTENOT_ASTOTA
1   2   INT_ASTENOT_ASTOTA
2   3   ER_ASTENOT_ASTTOM

It has more columns though.

Can you tell me how to read the column table_na line by line?

Example, reading by a txt:

with open(input()) as txtfile: #insert directory of txt
    x = txtfile.readlines()

But I have to read from this specific column only.

Also tried this (with cursor):

env.workspace = r"C:\Users\user\Desktop\05052\MyProj.gdb"

datasetList = arcpy.ListTables("*")

for dataset in datasetList:
     with arcpy.da.SearchCursor(dataset, "*") as cur:
          for row in cur:
              print row

But gives all that is in that table in one line.


texts:
I have a table (african_primates_17) set up in PostGIS and I have a view (finalised) generated from that. The view is then added into QGIS 3.
The simple idea was to overwrite the table by importing from an updated shapefile - and then expect the QGIS map to be miraculously updated also.
However, that did not work. I get the following error

Error 2
Creation of data source "public"."african_primates_17" failed:
ERROR:  cannot drop table african_primates_17 because other objects depend on it
DETAIL:  view finalised depends on table african_primates_17
HINT:  Use DROP ... CASCADE to drop the dependent objects too.
CONTEXT:  SQL statement "DROP TABLE IF EXISTS public.african_primates_17 RESTRICT"
PL/pgSQL function dropgeometrytable(character varying,character varying,character varying) line 15 at EXECUTE

Can anyone explain my error please? I thought the whole principle of views was for updates to be automatically made.


texts:
I want to use sharpmap to overlay a mapinfo tab file onto google or bing.

I have it working pretty much at the moment however the co-ordinates that bing/google show are in correct so my mapinfo file doesnt make any sense to teh region its being overlayed on. I think I need to transform the mapinfo co-ordinates to the bing/google type but I'm not entirely sure how to do this or if this is the correct thing.

Im using sharpmap 0.9

This is my code so far

myMap = new Map(new Size((int)imgMap.Width.Value, (int)imgMap.Height.Value));

var ogrProvider = new Ogr(@"C:\test.tab");

var layer = new VectorLayer("test");
layer.DataSource = ogrProvider;

var layerbing = new TileLayer(new BingTileSource(BingRequest.UrlBing, "", BingMapType.Aerial), "TileLayer");
//var layergoogle = new TileLayer(new GoogleTileSource(GoogleMapType.GoogleMap), "googlemaps");

myMap.Layers.Add(layerbing);
// myMap.Layers.Add(layergoogle);

myMap.Layers.Add(layer);
myMap.ZoomToBox(layer.Envelope);

When I run it up the "layer" has a SRID of 7 and the "layerbing" has a SRID of -1, I think this might be the issue and I need to transform the co-ordinates, but I'm not sure how to transfrom them and from/to what I need to.


texts: Matching GPS trackers POINT data between Nodes of OSM (Open Street Map)

texts:
I have GPS point data of vehicles from XYZ Tracking Company. The data has following columns:

UnitID
RecordDateTime
LATITUDE
LONGITUDE
Speed
Direction
Reason
Altitude

Data Dictionary is given below:

I want to aggregate (Median, Average) speed column with respect to  the given RecordDateTime between the nodes of OSM. So that i could get the Aggregated speed of all the vehicles on given Date and Time between the OSM (Open Street Map) Nodes. I want to prepare the above data like the example below?

So my summarized questions are:

. How can i match the lat/long between the NODES of OSM?

. And then how can i aggregate all the speeds on specific Date and time between the consecutive nodes of OSM in separate column like in example quoted above?

I want to do this through any of the following

GeoPandas,
Pandas,
PostGIS
PostgreSQL.

Looking for code or steps to perform this data preparation.


texts: greatest-n-per-group comparing geometry Postgis

texts:
I'm working on a query that I think could fit under the greatest-n-per-group category.

So here's a simplified version of the table:

So these are samples taken with a ~10 year span. So I have the exact same geometry stored in two or more rows.

I'd like to fetch only one row for each unique geometry, and the hard part is that I only want the latest taken samples (according to year column) and the data on that row.

I have tried something like this:

SELECT *
 FROM samples s
 INNER JOIN (
SELECT geom, MAX(year) as year
FROM samples
GROUP BY geom
) b ON ST_EQUALS(s.geom, b.geom) AND s.year = b.year

WHERE customer_id = 38330
ORDER BY provnr ASC;

This seems to give me a pretty good result, but it takes a lot of time. This specific customer has ~850 rows in a table of ~400k rows. That query takes about 7 seconds.
I have the geom and year columns indexed.

Is there a better way to write the query?
The query will be used in a real time app, so 7 seconds is way to long.

Really thankful for your'e help!


texts: How can I create an adapter between the ArcGIS REST interface and WMS or TMS?

texts:
I'm trying to use imagery available via the ArcGIS REST tile cache interface, however the tool I'm using, JOSM, only supports WMS or TMS. In addition the data from the ArcGIS server is in a "special" projection, Virginia Lambert Conformal Conic, which doesn't have an ESRI or EPSG wkid but does have WKT provided by the server.

Is there a way using free or open source tools to convert the tiles available from the server to WMS/TMS? I know this would involve downloading multiple tiles, then mosaicking and reprojecting them to create tiles suitable for WMS/TMS. I've come across MapServer and GeoServer, which seem very capable, but I wouldn't know where to begin, so any pointers would be helpful.

The server I'm interested in is here: http://gismaps.virginia.gov/arcgis2/rest/services/VBMP2006_2007/MapServer


texts:
I'm trying to filter Landsat images based on the percent of cloud cover for my ROI. I have tried to adapt code in link one  and link two but cannot get what I want. My ROI has been imported as a feaureCollection and covers 8 Landsat tiles.

This is my code

//select images from image collection
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1')
//LT05, LE07, LC08
.filterDate('1980-01-01', '2018-12-31')
.filter(ee.Filter.or(
ee.Filter.and(ee.Filter.eq('WRS_PATH', 208),
              ee.Filter.eq('WRS_ROW', 22)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 208),
              ee.Filter.eq('WRS_ROW', 23)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 208),
              ee.Filter.eq('WRS_ROW', 24)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 207),
              ee.Filter.eq('WRS_ROW', 22)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 207),
              ee.Filter.eq('WRS_ROW', 23)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 207),
              ee.Filter.eq('WRS_ROW', 24)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 206),
              ee.Filter.eq('WRS_ROW', 23)),
ee.Filter.and(ee.Filter.eq('WRS_PATH', 206),
              ee.Filter.eq('WRS_ROW', 24))));

var c = l8.map(function(img) { return img.clip(my_geom); });

var withCloudiness = c.map(function(image) {
var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
var cloudiness = cloud.reduceRegion({
reducer: 'mean',
geometry: my_geom,
scale: 30,
});

return image.set(cloudiness);
});
var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 20));
print(filteredCollection);

So I have two problems:

1: When using Landsat 8, the code runs, but it does not filter based on a threshold, 20% in this case. When i write the output to a file it gives two columns, cloud_cover and cloud_cover_land, both have observations of cloud greater than 20%.

2: When I try using Landsat 5 or 7 my file is empty...but I know there are scenes that are below 20% cloud cover for Landsat 5.

What I am I doing wrong?


texts: Satellite image vertically squeezed after being georeferenced

texts: I am using ArcGIS for Desktop 10.1. I am trying to georeference a Google Earth image in Christchurch area, New Zealand. I am using some CGPs I have got from this image (their Google Earth coordinates), but since the building I need to use to pick up the CGPs is still under construction it results only in the lastest satellite image and the resolution is not very good. The issue is that once georeferenced, the image is a bit squeezed vertically from the original one. I am using WGS84 as datum.


texts: QgsMapCanvas on-the-fly feature does not work

texts:
I am writing a simple program to display raster data using qgis api C++.
When I open raster data that has different coordinate systems, the on-the-fly feature does not work on my mapCanvas.
Can someone help me?

OpenGisData::OpenGisData(QWidget *parent)
: QMainWindow(parent)
{
    setupUi(this);
    QWidget* centralWidget = OpenGisDataClass::centralWidget;
    mMapCanvas = new QgsMapCanvas(centralWidget);
    QGridLayout* gridLayout = new QGridLayout;
    gridLayout->addWidget(mMapCanvas);
    centralWidget->setLayout(gridLayout);
    connect(actionOpen, &QAction::triggered, this, &OpenGisData::open);
}
void OpenGisData::open()
{
    QString fileName = QFileDialog::getOpenFileName(this, "Open data", "D:\\Test");
    if (!fileName.isEmpty())
    {
        QgsRasterLayer* rasterLayer = new QgsRasterLayer(fileName);
        if (rasterLayer && rasterLayer->isValid())
        {
            QList<QgsMapLayer*> layers = mMapCanvas->layers();
            layers.append(rasterLayer);

            mMapCanvas->setLayers(layers);
            mMapCanvas->zoomToFullExtent();
        }
    }
}


texts:
I have several points. These points are far from some lines.

I want to make the following:

If two points are closer to each other than the maximum distance of one of these points to the closest line, then connect the points to each other and connect the point that is closer to the line to the line. Is there a quick way of doing it?

I have used the xy to line to connect points to the line, but I do not know how to make this operation. Attached a snapshot to make it more clear.

(Edit - should the screenshot look more like this? The two points on the left are closer to each other than to the line, so they are connected to each other. The two points on the right are closer to the line than to each other, so they are connected to the line)


texts: QGIS 2.0.1 Dufour on Windows 8 64-bit: Processing toolbox always produces error saying output layers could not be opened

texts:
I have tried to use the Processing Toolbox in QGIS 2.0.1 a couple of times now (using SAGA and OTB) and in both cases the processing tools always return an error along these lines:

"Problem loading output layers: The following output layers could not be open(ed) - Output image: D:/Temp/output.tif

The above files could not be opened, which probably indicates that they were not correctly produced by the executed algorithm. Checking the log ... etc."

When I dismiss the error the original tool's window also closes, giving me no opportunity to check the log.

This is a brand new Windows 8 installation with a brand new QGIS installation (via the OSGeo4W online installer). D: is writeable and I have set up folders for the OTB toolbox.

What else should I do to get these toolboxes to work?


texts: Why does the Select by expression dialog not show field values?

texts:
I haven't found a similar question that already has an answer so I hope I can get one... Please forgive me if the answer is obvious.
I'm using QGIS version 2.14.0-Essen on my macbook pro with OSX Yosemite 10.10.5. I am trying the simplest select features using an expression

Here is the window I get:

I right-click the field (e.g. here "DIVISION") and select say Load all unique values but nothing happens. I can check layer values beforehand and enter the correct value and it works (e.g. "DIVISION"='CENTRAL') but it's cumbersome.

So my question is: How can I get this unique values window to show?


texts: Libraries for creating surfaces output in geotiff

texts:
I am looking for some libraries that can assist me in creating geotiffs from 3d polygons that becomes surfaces.

Given a set of polygons where each node is [x,y,z] coordinates, i would like to create a geotiff at a given resolution [r] that has elevation of the polygon at each sample point, such i can use it to subtract from DSM.

It is to be used to generate a surface that a DSM cannot penetrate. which can be found by subtracting the DSM from the geotiff surface generated from the 3d polygons above.

Do anyone have examples of libraries, it do not matter what the framework is - can be python or anything else.

It could also be combined to solve the real problem above such the intermediate step of creating a geotiff is not needed and it outputs the difference of the polygons to the DSM right away. Anyway, i need a way to represent 3d polygons and use them to subtract from DSMs.

(I dont want to be dependant on QGIS or other desktop installed software).

Can GDAL or PDAL solve the problem?


texts: Converting RGB image to data values using existing palette in FME

texts:
I have image (RGBA) data that needs to convert values using an existing palette. I have the palette available to get each RGB(A) -value translated to a correct value.
I need to convert RGB value into one band. After that I need to use palette and convert each RGB value to palette value.
Example I have red color pixel (RGB 255, 0,0) in original image. Then I have value 12 in palette which has correct value to convert red pixel to value 12.
How to do that?

The palette has total 25 values and looks like this:
RGB24
1 97,100,101
2 148,147,148
3 195,204,209
4 166,159,255
5 17,37,250
6 52,92,255
8 6,231,255
10 54,199,6
12 255,0,0
.
.

I have tried the followed FME process but its wrong.
Which transformers should be used in FME?


texts: Leaflet and tiler for R - can't get the correct tiles for displaying in Leaflet (wrong reference systeme or something)

texts:
I have got an issue with tiler (https://github.com/ropensci/tiler). I am setting up a webapplication and want to convert .grd files to tiles that can be read by Leaflet.
This example can be easily follow by using the files and putting them under a folder called /data, as you can see in the code, that is needed.
You can download the used files here at Gofile
I had to set the routes with tiler_options because there where some errors with systemvariables and the paths. But it runs and there is a output. But the tiles are obviously not in the right reference system.
This is my code (remember to install and load the package):
tiler_options(osgeo4w = "C:/OSGeo4W64/OSGeo4W.bat", python="C:/Users/____/AppData/Local/Programs/Python/Python37/python.exe")

tile(file= 'data/aoa_result_2.grd',
     tiles= 'layers/aoa_result',
     zoom="10-13",
     viewer= TRUE,
)

Even the passed viewer can't display the layer correctly. The XML file shows that the EPSG is 3857 but I guess that is not right for leaflet. It should be EPSG:4326 instead i guess.
<?xml version="1.0" encoding="utf-8"?>
    <TileMap version="1.0.0" tilemapservice="http://tms.osgeo.org/1.0.0">
      <Title>tmp_raster.tif</Title>
      <Abstract></Abstract>
      <SRS>EPSG:3857</SRS>
      <BoundingBox minx="7.40000709521776" miny="51.85465878258405" maxx="7.82132047624863" maxy="52.06491340494047"/>
      <Origin x="7.40000709521776" y="51.85465878258405"/>
      <TileFormat width="256" height="256" mime-type="image/png" extension="png"/>
      <TileSets profile="mercator">
        <TileSet href="10" units-per-pixel="152.87405654296876" order="10"/>
        <TileSet href="11" units-per-pixel="76.43702827148438" order="11"/>
        <TileSet href="12" units-per-pixel="38.21851413574219" order="12"/>
        <TileSet href="13" units-per-pixel="19.10925706787109" order="13"/>
      </TileSets>
    </TileMap>

I tried some reprojections but always got the same refsys.


texts: How does TileStache handle cache invalidation?

texts:
We are trying to set up our own openstreetmap server. From what I read here, you can use a variety of caches to speed up serving the data.

Now, once a week, the people of openstreetmap provides a differential update with the changes in the maps.

If I import those changes once a week, how does TileStache know that what it has generated is no longer valid? Does it handle it itself, or should I invalidate the cache from somewhere else?

Can I do it intelligently? ie, only invalidate the regions that contain changes.


texts: Query field disabled in edit data mode in pgadmin 4

texts: I want to make a SELECT statement in the "View/Edit Data" mode in pgAdmin 4 but the query field is disabled. How can I enable it? (I have all privileges for all, the database, schema and table and I am the owner of all as well)


texts: Calculate Distance Matrix for ALL possible combinations in the dataset in QGIS

texts:
I have a Vector Layer that is created out of a MySQL database and has the relevant points stored in WGS84 format. The import works well, however, when I try to calculate a distance matrix between all possible combinations of point via Vector->Analysis Tools->Distance Matrix I don't get my desired result. What I get with the "Linear standard matrix" option is a csv-File that contains all possible combinations for two arbitrary chosen fields with all other fields. When I choose "Standard distance matrix" Qgis crashes at 98% but the output file seems to contain the same arbitrary chossen fields in the first column (just multiple times) and all other fields in the first row.
What I would like to have is basically:

   ColumnID| Field 1 |Field2            |Field3
   Field1|-          |Dist Field1-Field2|Dist Field1-Field3
   Field2|-          | -                |Dist Field2-Field3
   Field3|           |                  |    -

Any idea how I do that? The format does not need to be in the same style but I need distances for all combinations.
Thanks


texts: Adding attributes from surrounding feature (exclave problem)?

texts:
Here is the problem:

I am creating polygons based on the result of closest facility tool in Arcmap. Polygons consist of municipalities with the same closest hospital (impedance is a travel time, not a spatial distance).

Since there are some municipalities with really close highway entrances, there are "islands" or "enclaves" in my polygons. I need these islands to be dissolved in surrounding polygons. For better understanding see the pictures below.

This is the example of a problematic municipality - "island" polygon

Result of the dissolve tool. Dissolve field from the input feature is Obec_spado.

My idea is to add a new field. If a polygon is completely surrounded by polygons with only one "Obec_spado" value (as seen in in the pictures 1/2; = "island" polygons), which is different from the "Obec_spado" of an "island" polygon, value of that new field would be "Obec_spado" from the surrounding polygons. In other cases value of the new field would be the same as original "Obec_spado".

Does anybody know how to do it? Or maybe any other solution? Since I need to do it many times with different input data and do not know Python, ModelBuilder advices would be great...


texts:
I'm trying to install the GEarthView plugin but without success. The following message pops up after installation:

"No module named zope.interface"

What is the problem?

QGIS 2.6.0
Windows7 (64Bit)


texts:
I have a simple script for processing thousands of rasters/shapefiles. i am not able to delete temporary files, which are created by processing. I can pass them to:

'memory:output_file'

but not for all of them, because sometimes I need to use gdal, which does not know memory. So I create gpkg/shapefile. But after 64 files (which are hold in memory, I get error from sqlite because it can not hold more than 64 open files...)

Sqlite ERROR: Too many connections: max 64

I am using standalone QgsApplication and calling processing.

The problem is, somehow the PyQGIS is locking files after processing.
Similar problem here:
QGIS: Release file lock on file used in a processing algorithm at the end of script

Another one here:
Releasing PyQGIS file locks?

Another one: https://issues.qgis.org/issues/12502 (which says fixed/implemented)

And many others. Their solution does not work for me. If i try to put it into:

QgsProject.instance().addMapLayer(layer)
QgsProject.instance().removeMapLayer(layer.id())
QgsProject.instance().mapLayers() ## which returns empty dict, all good here

And check there are no files in registry, there is still lock on files:
gpkg.wal and gpkg.shm

So i am not able to delete the file with:

os.remove(file)

My only solution is (which is not good, since I would have to call it all the time):

QgsApplication.exitQgs()

After that, I can delete the file. But if I try to start the application again, Python crashes.
Any ideas how to release the damn lock?

My script is eg:

def selectByLocation(input, predicate, intersect_lyr, output):
    process = "qgis:extractbylocation"
    params = {
        'INPUT' : input,
        'PREDICATE' : predicate,
        'INTERSECT' : intersect_lyr,
        'OUTPUT' : output
    }
    out = processing.run(process, params)
    return out["OUTPUT"]

def removeLowDepths(raster, threshold, rn, workspace = temp):
    workspace_memory = 'memory:'
    import uuid
    ## cant overwrite file because of lock, so i am using uuid to create many temp files, which i can not delete
    expression = "data < " + str(threshold)
    lowDepthMaskR =  setNull(raster, expression, os.path.join(workspace, "lowDepthMaskR.tif"))
    lowDepthMask = polygonize(lowDepthMaskR, workspace + f'//lowDepthMask_{str(uuid.uuid4())[:8]}.gpkg', epsg)
    fixed = fixGeom(lowDepthMask, workspace_memory +"lowDepthMask_fixed")
    lowDepthMaskSel = selectByLocation(fixed, 0, rn, workspace_memory + "lowDepthMaskSel")
    finalMask = selectByAttribute(lowDepthMaskSel, 'DN', 0, 1, workspace_memory + 'finalMask')

    return finalMask

Edit for comment:

def removeLowDepths(raster, threshold, rn, workspace = temp):
    workspace_memory = 'memory:'
    import uuid, time
    id = str(uuid.uuid4())[:8]
    expression = "data < " + str(threshold)
    lowDepthMaskR =  setNull(raster, expression, os.path.join(workspace, "lowDepthMaskR.tif"))
    lowDepthMask = polygonize(lowDepthMaskR, workspace + f'//lowDepthMask_{id}.gpkg', epsg)

    print(QgsProject.instance().mapLayers()) ## prints {} so empty dict
    del lowDepthMask ## se picture, still locked
    time.sleep(60)
    os.remove(workspace + f'//lowDepthMask_{id}.gpkg')
    return finalMask


texts:
I downloaded 10 .img files from garmin.openstreetmap.nl which are supposed to contain the topo tiles for the area I would like to add in qgis. I've been unable to find any instructions on how to use the files in QGIS.

Can someone explain or point me to where instructions are available to accomplish this?

I need to apply contour lines to a map of Thailand in qgis.
It appears the img tiles I currently have are unusable to achieve this.

Where can I download whatever file type is necessary containing the contour lines for Thailand and by what steps, using what application can that file be converted to one which can be added as a layer to qgis?


texts: Getting the Line or MultiLine coordinates in QGIS through PyQGIS

texts:
We can get point layer coordinates through the following code, but I want to know the coordinates of line or multiline. How to get these coordinates?
layer = QgsProject.instance().mapLayersByName("layerName")[0]
iter = layer.getFeatures()
for feature in iter:
  geom=feature.geometry()
  print(geom.asPoint().x(),geom.asPoint().y())


texts: Feather/buffer raster boundary during mosaic in QGIS

texts:
I have two rasters I would like to merge and also fill in no data cells along the boundary between the two rasters to create a smooth transition. Some of the cells along the boundary are no data cells, and they are not filled when I try the Mosaic raster layers tool by SAGA. I have tried both the feathering and blend boundary method for the "overlapping areas", but I suppose this only applies to the areas that actually overlap. Therefore, I'd like to feather/buffer my raster boundary to create data for the existing nodata cells for cells which are an x distance away from cells with data, similarly to how the built-in buffer tool for vectors works.
Is this possible with any built in tools in QGIS?
I've attached an image of my output from the Mosaic raster layers tool.


texts: I recently looked at some other posts relating to this topic but as far as I can tell there is currently not a working plugin or script for QGIS 2.0 to achieve a four color theorem for styling polygon vector data (it doesn't necessarily have to be 4 colors, just so that each neighboring polygon is a different color, though keeping it at a minimum number of 4-6 colors would be elegant).


texts: returning a LineString on selecting its midpoint

texts:
I have a  list of shapely LineString,

[<shapely.geometry.linestring.LineString object at 0x000001614D750CC0>, <shapely.geometry.linestring.LineString object at 0x000001614D750D68>, <shapely.geometry.linestring.LineString object at 0x000001614D750D30>]

and List of point geometry which are the mid-points of the Linestrings;

[<shapely.geometry.point.Point object at 0x000001614DC0D518>, <shapely.geometry.point.Point object at 0x000001614DC0D4E0>, <shapely.geometry.point.Point object at 0x000001614DCF9CC0>]

How to return Linestring on selecting its midpoint from the second list?


texts: Virtual Layers Clipping and Large Data Sets in QGIS

texts: When working with large vector data sets in QGIS is there a method where you can ask QGIS to only load a small subset of your data?  I know that this can be done based on attributes using a virtual layer, but I am wondering if it is possible to do this for multiple large datasets and to get the same geographic coverage.  Sort of like clipping the same geographic slice out of multiple data sets, except that I don't want to clip my data.  Is this possible.  I am not sure that I have the language to describe it.


texts: How to read coordinates in this specific format

texts:
I'm trying to create a shapefile consisting of lines from my python output. The output of the function is :

[[[400 421 427 114]]

[[227 475 237 356]]

[[408 383 427 124]]

[[456 322 500 327]]

[[436  68 441   8]]]

Where [X coord start point, y start point, x end point, y end point]

I've tried Shapely and Fiona but I can't get it to read this format.
How to iterate between the coordinates and putting them into a shapefile?


texts: Why is a single road sometimes split into separate layers?

texts:
I'm using data from the Transportation section of the National Map (US). They provide three layers for roads: Trans_RoadSegment, Trans_RoadSegment2, and Trans_RoadSegment3. I seem to recall this scheme from other GIS data--the same road is split into multiple layers.

In the example above, the road is split into three layers. Some of the segments overlap others, some take a slightly different route, and some do not overlap at all. In any case, if you were to display only one of the layers, there would be gaps in the road. Why is this done? Why isn't the road data all on a single layer? The National Map provides a data key poster, but I don't see any information about this split there.


texts: Obtaining "nClusters" value in Google Earth Engine

texts:
After training the kmeans model, how can i obtain the number of classes of the model? i've tried with get(), getInfo() and doesnt work.

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

Value that must to be extracted in "//var to_print_2 = clusterer.get('nClusters')"

// STEP 1: Import image we wish to classify
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1')

// Remove the clouds
var cloud_free = ee.Algorithms.Landsat.simpleComposite({
  collection: l8.filterDate('2018-01-01', '2018-12-31'),
  asFloat: true
})

// STEP 2: Draw your box which you will sample from!
// Make the training dataset.
var training = cloud_free.sample({
  region: geometry,
  scale: 30,
  numPixels: 5000
});

// Step 3: Create the clusterer and train it, play around with
// number of classes
var clusterer = ee.Clusterer.wekaKMeans(3).train(training);

print('clusterer',clusterer instanceof ee.ComputedObject);

var to_print = clusterer.getInfo();

print('clusterer_info',clusterer.toString());

print(typeof clusterer);

var clusterer_list = ee.List(clusterer);

print('clusterer_list',clusterer_list);

print(typeof clusterer_list);

//var to_print_2 = clusterer.get('nClusters');

//print(to_print_2);

//var to_print_3 = clusterer.Filter.filterMetadata('nClusters','equals',3);

//print(to_print_3);

// Step 4: Classify our image
var classified = cloud_free.cluster(clusterer);
print(classified);
var cluster_1 = classified.eq(0);
print('min',cluster_1);

// Display
Map.setCenter(-122.31040963183364, 37.782065864682096, 11)
Map.addLayer(classified.randomVisualizer(), {}, 'clusters');


texts: Converting vector cells to polygons in QGIS

texts:
I have a shapefile containing vector cells showing gradients of calcium pressure in soil. I want to convert the cell classifications into polygons in order to improve the visualization.
I'm new to GIS in general and even more so QGIS and have tried different things, like converting into points and then buffers, but without any satisfactory results.

Can anyone help?


texts: Why does clipping a raster alters cell size?

texts:
Using QGIS V 2.2.0 (OS-Windows 8), I need to prepare raster files (ESRI ArcInfo Grids) to  be later used in a species distribucion model using MaxEnt. Each of the original raster files (environmental variables for Guadalquivir, Spain) has a different spatial extension (but all share the same pixel size of 200, -200). However, to be able to use these layers in MaxEnt, all of them need to have the exact same spatial extent. Therefore, I clipped every raster using a polygon of Guadalquivir. After clipping, the new rasters have their original pixel size (200,-200) slightly altered (200.001,-199.942), which means trouble when running the maximum entropy algorithm.
The Layer Spatial Reference System for Guadalquivir is:
+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The Layer Spatial Reference System for the raster layers is:
+proj=utm +zone=30 +ellps=GRS80 +units=m +no_defs

I'm using the function raster >> Extraction >> Clipper followed by the options...

No data value = 0
Mask layer: Guadalquivir

Following, I save the new raster as GeoTIFF.
Questions:

Should the clipping be done using a polygon with the exact same projection? I am not sure to what extent projections differ in the layers above.

In case response to question 1 was "Yes" (same projections required), I did the clipping using both options: first by saving Guadalquivir as .shp with the raster's projection, and later the other way around, that is saving raster files with Guadalquivir's projection. Regardless of the option described in point 2), pixel size remains altered.

How can I keep the original pixel size of 200,-200?

I wouldn't bother too much about this seemingly slight difference if it weren't for the warning message that I get (see below) when trying to convert the GeoTIFF rasters to ASCII (a requirement for MaxEnt to work):
Warning 1: Producing a Golden Surfer style file with DX and DY instead
of CELLSIZE since the input pixels are non-square. Use the
FORCE_CELLSIZE=TRUE creation option to force use of DX for
even though this will be distorted. Most ASCII Grid readers
(ArcGIS included) do not support the DX and DY parameters.


texts: How do I specify an extent to my map without disabling zooming in?

texts: I'm a bit new to GIS in general so apologies for the newbish question but in ArcGIS 10 I'm trying to get it so my map is only for my local city government (where all my data is) and a small surrounding area. I'd like to have it set up so that by dragging the scroll bars I'm confined to that area instead of the current scrolling around the entire earth. It renders the scroll bars useless besides a quick arrow nudge here or there. I've tried different extent settings on my data frame from this page from ESRI to no avail. I want to still be able to zoom in to inspect features closer but all viewing extent settings seem to disable zooming in and out. Is my base map causing a problem since it's worldwide imagery? All my other data is confined to my city and some small surrounding area. Am I out of luck in doing this? Thanks in advance!


texts: Proxy Error while using Google Earth Engine Python API

texts:
I am trying to use Google Earth Engine (GEE) Python API using the following code:
import ee

#ee.Authenticate()

ee.Initialize()

With the above code I am successfully able to Authenticate myself, but in Intialization step I am getting the following error:
ProxyError: ('Cannot connect to proxy.', OSError('Tunnel connection failed: 407 Proxy Authentication Required'))

I am getting the above error as I am trying to access this API from behind the firewall system. I have added proxy details ("http_proxy", "https_proxy") in my system environment variable in this format 'http://username:password@proxy:port' but still I am getting the above error.
Can someone please suggest me how to overcome this above error.


texts:
I have a MapInfo 11 user who wishes to use ArcGIS 10.3.1 colour ramps. I'd like a seamless method to get from A to B.

So far I've screen grabbed the ramp and saved it as a raster, exported it to points,  calculated the raster's RGB values into the attribute table, and run a frequency on the table.. This gave me a list of all the RGB combinations in the ramp. Apparently MapInfo users can import RGB values.


texts: Downloading Google Earth images in universal map downloader

texts:
How can I download Google Earth images in universal map downloader in 2015 and 2020?
There isn't any specific option for select year.


texts:
I am working on a personal web gis client and the first thing I am implementing is support for OGC WMS 1.3.0.

Right now, I am in the middle of greying out layers in my layer tree which are not visible at the current view scale.

From the official documentation I read:

  The  and  elements define the range of scales for which it is
  appropriate to generate a map of a Layer. Because maps from a WMS may be viewed on arbitrary displays rather
  than on paper, the values used are actually the scale denominators relative to a common display pixel size. (...) For the purposes of this International Standard, the common pixel size is defined to be 0,28 mm × 0,28 mm.
  Because arbitrary clients can request maps from a server, the true pixel size of the final rendering device is
  unknown to the server.

Probably is either that English is not my first language and that I am very bad with numbers, but I am struggling to figure out the meaning of this sentence.

Take this WMS for example. Some layers have <MaxScaleDenominator> equals to 94494.047619.

First of all, how can I read this number in terms of scale? Does it represent the common 1:94494.047619 scale, or do I have to do some calculations with the standardize 0,28 mm × 0,28 mm pixel size to understand it?

Also, if

  the true pixel size of the final rendering device is unknown is there
  a way to to the server

is there a way to retrieve it client side (if that is indeed required for what I need to do)?

EDIT1

So I already got pretty good results. I have the below code which logs to console the current scale of my view (might be not 100% clean, but it seems to work, I adapted it from here):

import {METERS_PER_UNIT} from 'ol/proj/Units.js';

// return current scale in map units according to DPI
function mapScale(map, dpi) {
    var unit = map.getView().getProjection().getUnits();
    var resolution = map.getView().getResolution();
    var inchesPerMetre = 39.37;
    return resolution * METERS_PER_UNIT[unit] * inchesPerMetre * dpi;
}

export function getScale(map) {
    map.getView().on('propertychange', function(e) {
        switch (e.key) {
            case 'resolution':
                // log DPI to console
                console.log("DPI:" + window.devicePixelRatio);
                // log scale in map units to console
                console.log("1:" + mapScale(map, window.devicePixelRatio * 96));
                break;
        }
    });
}

Now I can see in the console when I am farther than the aforementioned <MaxScaledenominator> (94494.047619), and indeed, if I turn on the layer it won't show up till I have reach at least this value (actually it turns on at 1:72223 at minimum as I am using OpenStreetMap basemaps and this is one of its fixed scale level).

This seems to confirm that the numbers are meant like "traditional" scales (i.e. 1:94494.047619). Can anybody please confirm and elaborate a bit on this?


texts: QGIS 2.18 Las Palmas Error Screen upon project opening. Win 10

texts:
I am working on a trail map for a Non-Profit.  Created a few maps last couple of years with previous versions of QGIS and just updated to 2.18.  Running a Windows 10 OS.

QGIS loads fine, but when project starts to load get pop up error message with cryptic login requirement with unknown user and password already in place.  Of course this is Not valid and will not clear.  I do not know for sure what site it is trying to log in to.  Program will not load and cycles back to error message.

To the best of my ability it appeared to be trying to access the qgiscloud.com.  I did change the default password being displayed with my qgis cloud password to no avail.  Have attached image of error message.


texts: QGIS Print Composer does not change from landscape to portrait orientation?

texts:
Using 2.10.1, I open a new Print Composer and immediately change the Page Setup Orientation from Landscape to Portrait.  But the orientation of the screen does not change; it continues to display landscape.

Likewise, changing Paper Size does not affect the size/orientation of the Print Composer window.

If I add a map to the Print Composer and Export as PDF (regardless of what Page Setup is set to), the result is always landscape.

It appears that Print Composer is only capable of A4 landscape.

What am I missing?


texts: How to use old WKT(version1) coordinate reference system in ArcGIS Desktop?

texts:
I have a local coordinate reference system defined with WKT:

            COMPD_CS[
                "test",
                FITTED_CS[
                    "",
                    PARAM_MT[
                        "Local_Grid",
                        PARAMETER["E0",131473.823],
                        PARAMETER["N0",25419.216],
                        PARAMETER["K",0.999997513402],
                        PARAMETER["DE",-522111.156],
                        PARAMETER["DN",-2500318.967],
                        PARAMETER["Beta",359.022486769444]
                    ],
                    PROJCS[
                        "",
                        GEOGCS[
                            "",
                            DATUM["WGS 84~1",SPHEROID["xian80",6378137,298.3] ,TOWGS84[0, 0, 0, 0, 0, 0, 0]],
                            PRIMEM["Greenwich",0],
                            UNIT["Radians",1],
                            AXIS["lng",East],
                            AXIS["lat",North]
                        ],
                        PROJECTION["Transverse_Mercator"],
                        PARAMETER["latitude_of_origin",0],
                        PARAMETER["central_meridian",114],
                        PARAMETER["scale_factor",1],
                        PARAMETER["false_easting",500000],
                        PARAMETER["false_northing",0],
                        UNIT["Meters",1],
                        AXIS["east",East],
                        AXIS["north",North]
                    ]
                ],
                VERT_CS[
                    "",
                    VERT_DATUM["Ellipsoid",2002],UNIT["Meters",1],AXIS["Height",Up]
                ]
            ]

My first try is to save it as .prj file and open with ArcGIS, but ArcGIS cannot recognize it.

After reading the official doc (http://www.geoapi.org/3.0/javadoc/org/opengis/referencing/doc-files/WKT.html , http://docs.opengeospatial.org/is/12-063r5/12-063r5.html), I knew it defines a transformed projection coordinate system with old WKT standard (version1), and some concepts in it is no longer used in current version ("COMPD_CS","FITTED_CS").

My second try is to search how the transformation works (PARAM_MT[..] part) ,hoping to rewrite it in current WKT standard or make a transformation in ArcGIS, but found no useful information.

Does anyone know how to make use of this old WKT defined CRS?


texts: Shapefile created in QGIS does not fit in eCognition project

texts:
I have various shapefiles - some created by others, some created by my self. They should all have the same size of the underlying Raster file.

The shapefiles have been created on top of a Raster file (in QGIS) with Coordinate System WGS84 /UTM zone 45N. I checked the Options in QGIS - Shapes and Raster use the same Coordinate System.

When importing the Shapefiles to Erdas Imagine, everything looks fine.

In eCognition the shapefiles doesn't fit to the extend of the Raster file.

Where can I solve this problem? Do I have to change something in QGIS or is it a eCognition problem?


texts: How to create custom coordinate system in ArcGIS Desktop?

texts:
I have a series of control points with xy coordinates in both NAD83 UTM and in a local grid.

How can I go about creating a new projection file in ArcGIS 10 for the local grid?