texts: Adding customized basemap ArcGIS vector tile to QGIS

texts:
I customized a basemap through the platform ArcGIS Vector Tile Style Editor (beta) and I would like to share it and use it also with other software like qgis.

I am using this link:
https://www.arcgis.com/sharing/rest/content/items/15a61dd0045f40b0836dc73af2ed9d3d/resources/styles/root.json

Does somebody know if customized vector tile created in arcgis can be usually added to QGIS or other software?

What do I need to do?


texts:
Using ArcMap 10.3, I am trying to perform map algebra on a raster with 2013 annual PRISM precipitation for the U.S. (http://prism.oregonstate.edu/recent/). The raster has 1 band an is floating point.

Using Spatial Analyst Tools > Map Algebra > Raster Calculator (i.e. not within model builder), I have tried to create a new raster by executing these expressions (separately):

"PRISM_ppt_stable_4kmM3_2013_asc.asc" * 0.1 and Float("PRISM_ppt_stable_4kmM3_2013_asc.asc" * 0.1)

Both expressions yield this error:

Raster(r"PRISM_ppt_stable_4kmM3_2013_asc.asc") * 0.1
ERROR 000539: Error running expression: rcexec()
Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 5, in rcexec
  File "c:\program files\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 4333, in Times
    in_raster_or_constant2)
  File "c:\program files\arcgis\desktop10.3\arcpy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "c:\program files\arcgis\desktop10.3\arcpy\arcpy\sa\Functions.py", line 4330, in Wrapper
    ["Times", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 999998: Unexpected Error.

Failed to execute (RasterCalculator).

What am I missing here?


texts: Can't install QGIS 2.6.0 need SQlite3 on OSX Yosemite

texts:
I need to install QGIS on my mac with 10.10 Yosemite, I have already installed:
- SQLite3_Framework-3.8.7.3
- GDAL_Framework-1.11.1-2
- GSL_Framework-1.16-1
- NumPy-1.8.0-1

but when I launch QGIS.pkg it says it requires SQlite3 framework. Please help me


texts: Can't connect with my PostgreSQL database with ArcGIS

texts:
I have created a database in PostgreSQL, spatially enabled it, and now I need to connect it with ArcGIS (we use only ArcGIS products here on my university), everything on localhost, I've tried everything I could find here on stackexchange and on the ESRI support forums, nothing seems to work, I'm a little desperate

Here, I will describe my entire situation and the solutions I've tried:

I used to have a 64bit PG instalation on my 64bit machine, reinstalled both Postgre and PostGIS and recreated a database on 32 bit according to ESRI support forums, nothing changed.
I use 10.5 ArcGIS Desktop, so, it comes with the driver (client files) for Postgre with it according to another post here in stackexchange, didn't work, I've also tried to replace those 5 .dll's files from PostgreSQL folder in the ArcGIS\bin folder, nothing changed.
When I had a 64bit installation I have also tried to copy the "C:\Program Files (x86)\ArcGIS\Desktop10.5\DatabaseSupport\PostgreSQL\9.5\Windows64" files to my PostgreSQL folder, when 32bit I did not, first because the last version there is 9.5 and my PG is 10.4 and then because it seemed like those files were for 64bit (the Windows64 at the folder's name)
I have changed the pg_hba.conf to accept any connection as above

IPv4 local connections:
host    all             all             0.0.0.0/0            md5
The postgresql.conf file is set to listen to any connection (at least I think so), configuration is set as follows

listen_addresses = '*'
Just to note, anytime I did any alteration on the conf files I had the postgresql service restarted
I can connect with pgAdmin and QGIS locally and remotely JUST FINE (frustating how fine it is)

Here I'll have you some prints of my sad situation

That's how my pgAdmin looks when showing the properties of my connection

Now, on ArcCatalog, when I try Database Connection and try to open the droplist of databases

And if I'm stubborn enough and try to manually set the database name

And then, it won't connect either.


texts:
I have created a gxp client where I have added some layer in the initial map settings and some layer can be added using gxp Add Layers plugin. When I add these layers to the map there is no issue with them, I can get their properties in gxp layer properties plugin, print and getfeatureinfo for them. But I also need to add some of the layers using and apply filter on it. Now the problem in this method is that it doesn't create GeoExt.data.LayerRecord and these funtions are not available for this layer. How can a create a layer record for the following layer so that getfeatureinfo, print and layer properties function should work.

filterLayer = new OpenLayers.Layer.WMS(title,
        "/geoserver/wms", {
        layers : 'ws:layername',
        transparent : true,
        styles : style
    }, {
        isBaseLayer : false
    });
app.mapPanel.map.addLayer(filterLayer);

I also tried but this also didn't work.

var recordType = GeoExt.data.LayerRecord.create();
var record = new recordType({layer: filterLayer, title: "test Layer"});
app.mapPanel.layers.add(record);


texts: Maxent output clips to minimum extent of all input variables

texts:
I am experiencing an issue with the spatial extent of my output ASC files. I am using the UK Land Cover Map (LCM) 2007 at 25 metres resolution as a categorical variable in my Maxent modelling, alongside select continuous Bioclim variables interpolated (spline) to 25 metres and to the same extent. However, as the Bioclim files are originally at 1 km resolution, in some areas along the coastline they do not fully cover the sand dunes and mudflats covered the LCM which I am trying to model for.

The 'extent' of the Bioclim layer is the same as for LCM but these missing cells are given a 'No Data' value in the input ASCs. The output Maxent ASC extents are clipped to that of the Bioclim layers but I want them to extend to the extent of the LCM (which is included in the analysis and is in fact the 'most important' variable). Coastal sand dunes are an important habitat for my species and omission of these areas will negatively affect the accuracy of further connectivity analyses based on my Maxent outputs.

Is there a setting in Maxent to resolve this, or would I need to interpolate the input Bioclim layers outside of their original extent? If so, would this method have merit (bearing in mind that Bioclim is already an interpolated surface)?


texts: Define or convert projection when netcdf saved as tif raster?

texts:
A netcdf data structure showed as follows:

Dimensions:  (lat: 360, lon: 720)
Coordinates:
  * lon      (lon) float32 -179.75 -179.25 -178.75 -178.25 -177.75 -177.25 ...
  * lat      (lat) float32 -89.75 -89.25 -88.75 -88.25 -87.75 -87.25 -86.75 ...
    time     datetime64[ns] 1901-01-16
 Data variables:
    cld      (lat, lon) float64 nan nan nan nan nan nan nan nan nan nan nan ...
    stn      (lat, lon) int32 -999 -999 -999 -999 -999 -999 -999 -999 -999 ...

I want to save this layer as .tif raster with wgs84 proj. How to convert the projection?


texts: Not matching schema problems when appending list of CSV's into existing file geodatabase table?

texts:
So far I have the following code for creating the FGDB table and defining fields which works just fine. But when I try to append a list of csv's to this table I always get a matching schema error.

Do you have any idea, what may be wrong? Code and csv structure posted below.

import arcpy
import os

arcpy.env.workspace = "C:\\Users\\SAMSUNG\\Desktop\\python_script"
gdb = "FileDatabase.gdb"

if arcpy.Exists(gdb):
    print "Database exists."
else:
    print "Database doesn't exist."
    fgdbPath = "C:\\Users\\SAMSUNG\\Desktop\\python_script"
    fgdbName = "FileDatabase.gdb"
    arcpy.CreateFileGDB_management(fgdbPath, fgdbName)
    print "Database has been created."

tabPath = "C:\\Users\\SAMSUNG\\Desktop\\python_script\\FileDatabase.gdb"
tab = "Table"
fullPath = tabPath + "\\" + tab

if arcpy.Exists(fullPath):
    print "FGDB table exists."
else:
    print "FGDB table doesn't exist. Table will be created."
    arcpy.CreateTable_management(tabPath, tab)
    arcpy.AddField_management (fullPath, "transect", "LONG", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "point", "LONG", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "Xcoord", "DOUBLE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "Ycoord", "DOUBLE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "date", "DATE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "depth", "LONG", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "paramA", "DOUBLE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "paramB", "DOUBLE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "paramC", "DOUBLE", "", "", "", "", "", "REQUIRED")
    arcpy.AddField_management (fullPath, "notice", "TEXT", "", "", "", "", "", "REQUIRED")
    print "Table has been created."

For appending listed csv's I'm using the following script:

import arcpy

arcpy.env.workspace = "C:\\Users\\SAMSUNG\\Desktop\\python_script"

csvList = arcpy.ListFiles("*.csv")
arcpy.Append_management(csvList, workspace + "\\" + "FileDatabase.gdb" + "\\" + "Table", "TEST")

The csv files are without header and look like this:

1,5,504577.885850,1158507.476675,15.04.2015,40,3.56,0.007,23.34,"text"
1,6,504577.885850,1158507.476675,15.04.2015,40,3.56,0.007,23.34,"text"


texts: Is there a way to prevent labels from being on top of everything else?

texts:
I'm having a problem in QGIS 1.7.4. The labels I have produced for one layer are overlapping an entirely different layer that has a higher precedence. This should not happen. Am I doing something wrong or is this a bug?

And just to make sure my question was clear, I'm asking if there is a way to prevent labels from being on top of everything else regardless of layer order


texts: Connecting to ArcGIS ImageServer in QGIS?

texts:
I'm having trouble opening an ArcGIS ImageServer layer in QGIS. The city of Fort Worth, TX has a 2017 ortho imagery service which can be opened in ArcGIS. I've tried to connect to it as a WMS Layer and as an ArcGIS MapServer Layer through QGIS, both of which prompt a bad request error message.

This is the link that shows in the ArcGIS layer properties:

https://mapit.fortworthtexas.gov/ags/rest/services/Imagery/Orthos_2017/ImageServer


texts:
I try to add GeoJSON layer to my openlayers 3 by the following code.

        new ol.layer.Vector({
            title: 'added Layer',
            source: new ol.source.GeoJSON({
               projection : 'EPSG:3857',
               url: 'poly1.json'
            })
        })

Getting the following error,
Uncaught TypeError: ol.source.GeoJSON is not a constructor


texts: How can I make this model more efficient?

texts: The issue I am running into is when the amount of data scales up the model takes a very long time to process.  The step that seems to be bogging down the model is the extract by mask portions while the rest of the model seems to run at about what I would expect.  My thought is that there is some processing extents or some environmental setting I am missing.  A snip of the execution shows the rough times it takes to process each step.


texts: Creating yyyy-mm-dd column in attribute table from Julian day in QGIS?

texts: I am using the time manager for a set of data and I need to create a column in my attribute table of the date in the format yyyy-mm-dd.  I already have a column in my attribute table that contains Julian days (e.g. values of 1, 2, 3, etc. that I want to correspond to 2010-01-01, 2010-01-02, 2010-01-3, etc. respectively).  I've tried using the field calculator to do this using Converting from Julian Date to Calendar Date? but I was unable to get an output.


texts: How to query OpenStreetMap data served by GeoServer?

texts:
We have a specific point and would like to count all features in a circular area with center at that specific point.

The data is:

Stored in a shapefile layer created in QGIS
Served by GeoServer
Browsed by a web browser using OpenLayers


texts: Using arcpy.CreateMapTilePackage in ArcGIS Pro?

texts:
I have the  arcpy.CreateMapTilePackage_management() script working however I want to have the extent of which the tool tiles set to a layer and I'm having a bit of trouble understanding the syntax. Below is a link to the ESRI documentation:

http://pro.arcgis.com/en/pro-app/tool-reference/data-management/create-map-tile-package.htm

Here is a small bit of code with example code I wrote. I'm hoping someone can chime in and guide me in the correct direction on how to use the correct syntax on getting the tool to use the extent of a layer.

arcpy.CreateMapTilePackage_management(mxd_file_name, "ONLINE", lyr.datasetName + ".tpk", "PNG8", "10", EXTENT = <Conductor>)


texts:
I have just installed QGIS 1.8 on Windows 7.  Once the install completed and I tried to open QGIS (via the desktop icon) nothing happened.  I had this issue when installing QGIS 1.7.  I was able to get 1.7 to open by using the qgis.bat file, but it is not in the the bin folder in the QGIS 1.8 version.  I began to have this issue after installing OpenGeoSuite 4W.  I thought I uninstalled everything I needed.  I am worried I may need to re-install Windows.  Does anyone have any ideas as to what I can do to get QGIS 1.8 to load and open.

Thanks...


texts: Georeference Raster File with WorldFile into an aerial photo using GDAL

texts:
Possible Duplicate:
  How can I transform a georeferenced JPG to a georeferenced TIF?

UPDATE:Corrected link

I must note that I am new to GIS so pleas be understanding.

PROBLEM:

I have a raster file in this case it is an ortomap:
http://imageshack.us/photo/my-images/703/ortofoto1.jpg/

and it has JGW (World File)

0.683772574816354
0.000000000000000
0.000000000000000
-0.683739485024190
395057.114562241190000
741149.458682749890000

Now I want to use GDAL to apply the jpg+jgw to this aerial image:

http://imageshack.us/photo/my-images/688/aerialphoto.jpg/

so that it can have the coordinates according to the jgw file of the jpg ortomap.

What gdal console commands should I commit on those 3 files and what parameters should I consider?
(I was thinking gdal_translate/gdalwarp maybe)


texts: How to add logical operators for WFS in OpenLayers2?

texts:
How I can add logical operators AND, OR both together for WFS in OpenLayers.

When I try to run above code in browser nothing is being displayed and Firebug is showing an ReferenceError: init is not defined init() error.


texts: How to get closest Vertex value from Selected Feature by point Using qgis python

texts:
For example
I have a point ( 0 , 0 )
And I have a selectedFeature
How can I get value ( x , y ) of closest vertex from selected Feature by Using Python Console in Qgis

image for example :


texts: Reclassify Tool Resulting in NoData where Input has Data

texts:
Software: ArcGIS Pro 2

I am having an issue with the Reclassify tool. I am attempting to reclassify raster data into three classes so that I can do a habitat suitability analysis. When I use the Reclassify tool on some raster data I have (not all raster data I have) the resulting reclassified raster has several NoData strips in seemingly random areas across my study area. These strips are not NoData areas within my Input raster. Please see examples below:

Input Raster:

Reclassified Raster (white = NoData):

I haven't been able to figure out exactly what could cause this, mainly because I'm not sure where to start. Has anyone experienced this issue?


texts: Searching Geodatabase for Feature Class or Table by its name?

texts:
I'm looking for a way to search Geodatabases for a Feature Class or Table by using its name.

Is there any function like this available?

Obviously windows search only returns the random ID's for the tables within a Geodatabase. I'd like to actually search for datasets within them.

I think this would help a lot of people, because I know I'm not the only one out there who stores spatial data in Geodatabases, only to archive it with a thousand other datasets and not want to sift through it.

I'd like to search these Geodatabases within ArcCatalog.


texts: Timemanager not working with Integer Timestamp format

texts:
I'm trying to create an animation using Timemanager on QGIS 2.18 but I'm struggling.

I have a dataset composed by a list of businesses and their opening/closing time. This is what the dataset looks like:

'Start' and 'End' columns contain the day and the opening hours of the business converted in Unix Timestamp.
I upload my CSV in QGIS and add the layer on TimeManager. The plugin seems to accept the data. This is the slider after I uploaded my data:

Even if everything looks ok my animation doesn't work, nothing happens to the points on my map. I simply need the points to appear when the business is open and disappear when the business is closed.
Do I do something wrong? I never used this plugin before, so I don't know how to tackle the problem.


texts: pgrouting: how to temporarily increase some edges' cost without affecting concurrent queries?

texts:
How to temporarily increase some edges' cost for one query without affecting other concurrent query?

The case is, I would like to develop a system that allows users to query the shortest path of a road network. Users are able to avoid some roads, what I think is to increase those road's cost temporarily for rerouting. But if I update the table directly, it will affect the original data and other concurrent queries.

Should I just use transaction to isolate the data? Is there any better way to do so?


texts: Raw State Data on Health and Social Outcomes

texts:
I am looking for raw health and social outcome data (poverty, obesity, education levels, etc.) by state. I would like to download it as a csv or excel file.  It needs to have the state fips listed in one of the columns.

There is a lot of county data out there, but can anyone point me in the right directions for state data?


texts:
I am attempting a simple data extraction from an Oracle database to a CSV file using ogr2ogr.exe
(the version packaged with QGIS 3.8 i.e. ..QGIS_3.8\bin\ogr2ogr.exe )
I have a field ("Mapsheet") with values like : '092H034','114P055','082E033'
Values with the letter 'E', bracketed with integers (ex. '082E033') are a problem - ogr2ogr extracts these as exponential notation! What I see in my output CSV file is a float, not a string: 8.2e34
I have employed the following SQL as a workaround:
select (s.mapsheet || '_') as Mapsheet from myTable as s
That SQL is contained in a separate myFile.sql file, so the call looks something like this:
ogr2ogr.exe -f "CSV" test.csv OCI:user/password@DBP01 -sql @T:\myFile.sql -overwrite
I would of course prefer to cast the data using OGR SQL, since my workaround above appends an extra '_' character to the end of the value ... not ideal.
I've tried these CAST statements but they still result in a float number (ex. 8.2e34) being written to my CSV file.
cast(s.mapsheet AS char(7)) as Mapsheet
cast(s.mapsheet AS character(7)) as Mapsheet
So, how can I properly cast this field to be string values?


texts: Monthly Water Area Monitoring in Google Earth Engine?

texts:
I've computed water area for each 8 days for about 20 years. I don't know how can I compute monthly area for these years? or How can I change 8 days product to the monthly in this code?

code link: https://code.earthengine.google.com/590d4cf09f90e591d69ed3a03039d4f7

Map.centerObject(geometry)

var modis = ee.ImageCollection('MODIS/006/MOD09A1')
.filterBounds(geometry)
.filterDate('2000-01-01','2020-01-01')
.select('sur_refl_b03','sur_refl_b02')
.map(function(img){

  var id = img.id();

  var scale = img.multiply(0.0001);
  var clip = scale.clip(geometry);
  var index = clip.normalizedDifference(['sur_refl_b03','sur_refl_b02']);

  var thr = index.gte(0.1);

  var mask = index.updateMask(thr);

 var area = mask.multiply(ee.Image.pixelArea().divide(1000000.0));

  return area
  .copyProperties(img,['system:time_start','system:time_end'])
});

print(modis);

var chart = ui.Chart.image.series(modis, geometry,
ee.Reducer.sum(), 1000, 'system:time_start');

print(chart);


texts:
I have a code that analyze NDVI for any shapfile and then compute the mean and the standard deviation for each image.
My goal is to be able to do more calculation with the statistics I have calculated, but seems like because it's a dictionary I can't calculte.

I want to calculate the next things-

Standard deviation/2
mean+ (standard deviation/2)
mean-(standard deviation/2)

I have tried to calculate it by I get ''NaN'' all the time.

this is the code I have written-

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2019-06-01', '2019-06-30')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B1','B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

//test if clipping the image collection worked
Map.centerObject(geometry,9);
Map.addLayer(clippedCol.median(), rgbVis, 'RGB');

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

//analyze images from image collection collection
var listOfImages = withNDVI.toList(withNDVI.size());
var listOfNumbers =[0];

for (var i in listOfNumbers) {
  var image = ee.Image(listOfImages.get(listOfNumbers[i]));
  var meanDictionary = image.reduceRegion({
  reducer: ee.Reducer.mean(),
  geometry:geometry.geometry()
  });

 var STDDictionary = image.reduceRegion({
  reducer: ee.Reducer.stdDev(),
  geometry:geometry.geometry()

});

print(i,'mean',meanDictionary,'standard deviation',STDDictionary);

 }

var std2=STDDictionary/2;

My end goal is to be able to classify the image into 3 classes-

pixels with value that is between mean-std/2 and mean+std/1
pixels with higher value than mean+std2
pixels with value that is smaller than mean-std2


texts: Are there any embedded spatial databases for Java? I've come across H2 Spatial, but it was difficult to tell the status of that and if it can be used as an embedded spatial database.


texts: gdal_translate- why jpg compressed tif is 2 times greater than jpg file?

texts:
Using the following two commands yields completely different file sizes:

gdal_translate test.tif -of GTiff -co "COMPRESS=JPEG" -co "JPEG_QUALITY=100" output.tif

and

gdal_translate test.tif -of JPEG -co "QUALITY=100" -co "WORLDFILE=YES" output.jpg.

(I Wrote the commands from my memory, not sure that the syntax is correct)

The jpg file (second command) is 2 times smaller than the compressed tif file. Why?


texts: How to merge rows in the attribute table while adding up one column

texts: I intersected two layers and now I would like to merge the rows that have the same id while adding up the last column using QGIS.


texts:
Python newbie needing help with a simple task.  I'm trying to create a list of stand-alone tables in an sde workspace.

Code below works just fine to list feature class names and paths and table names.

fcList = gp.ListFeatureClasses()
for fc in fcList:
    descFC = gp.Describe(fc)
    print "FC Name: " + fc + ", Path: " + descFC.CatalogPath

tblList = gp.ListTables()
for table in tblList:
    print table

However, I am not able to use the Describe method to grab the tables' catalog path.

tblList = gp.ListTables()
for table in tblList:
    print table
    descTBL = gp.Describe(table)
    print descTBL.CatalogPath

I think my issue is in the gp.Describe line.  What's wrong with gp.Describe(table)?  From looking at code samples fc = feature class and fd = feature dataset.  What abbreviation is used for table?


texts: Transforming coordinates between wgs84 and UTM ED50 31N

texts:
I must tranform a point cloud from wgs84 (geographic) to UTM (projected) ED50 for zone 31North.

I have several free webs who offer me this service, and the answer is slightly different on every one of them.

Worse of all, I transformed it using ArcCatalog 9.2 ToolBox (the one you find on 'Data management tools->Projections and transformations->Feature->Project'.

Finally I found a little excel tool on this link http://www.aristademidi.es/app/download/5793157100/Conversi%C3%B3n+GPS_UTM.xls
(i think it can be credited to geographer Gabriel Ortiz from www.gabrielortiz.com , but this detail could be inexact) and quite to my despair, it turned out that the locations this tool gives seem more exact (when I draw it over my cartographic base) than the ones from Arcgis. The difference between this and the excel tool is systematically a bias +10m E.

On ArcCatalog I used the parameters from WGS 1984.prj to European
Datum 1950 UTM Zone 31N.prj, using something called transformation ED_1950_To_WGS_1984_31.
On the excel tool I used the upper side of the second sheet.

Where is absolute truth on this matter and where could I find the 'real' numeric transformations in case I must code a tool that I can really trust?


texts: How to install CartoDB driver for ogr2ogr?

texts:
I'm trying to import data from my postgis DB directly into CartoDB. I've seen in this forum's post that I can do that with ogr2ogr, a GDAL tool, in this way:

ogr2ogr -t_srs EPSG:4326 --config CARTODB_API_KEY XXXX -f CartoDB
"CartoDB:account" PG:"host=localhost user=postgres
dbname=cartodb_dev_user port=5432" "public.table_name"

However, when I try, I get this error message:

Unable to find driver `CartoDB'.

I don't have the driver for CartoDB installed, but I'm not sure how to install it. I see there's this help site, but it doesn't mention how to install it.

My current version of GDAL is 1.10.1 and I'm working on Ubuntu 14.04 machine.


texts: Could $area in field calculator be in cm?

texts:
I know there are several questions regarding issues with with $area function in the field calculator out there, however, but I could not find anything on this exact topic.

All my layers and the projects are set with CRS NAD83 UTM Zone 7. I recently downloaded the new version (2.18) of QGIS.

My $area output for the layer is huge (2241602779.0000). I know the outputs are supposed to be in meters, however, if they were in centimeters it may make sense.
Is this possible? Could I have some setting wrong in my field calculator?

I tried starting a new project, setting project projection and only loading the layer of interest with the same results.

I looked at the units for distance and area measurements in the project properties, and these are both set in meters, and square meters.

I am not certain how to proceed.


texts: Does QGIS have Generalization toolset like ArcGIS Desktop?

texts:
In ArcGIS Desktop there are Generalization tools like:

Merge Divided Roads
Collapse Road Detail
Resolve Road Conflicts

Does QGIS have this toolset or is it possible to run these algorithms on QGIS?


texts: Get nearest line to point with v.distance in QGIS

texts:
I want to retrieve the nearest line to a point with the GRASS function v.distance(). But I think I do not fully understand how to use this function. What I have are two layers (in the same CRS -> EPSG:31255): one point-layer which contains the locations of all my sensors and one line-layer which contains a streetnet. The following picture shows both my layers:

With v.distance, I want to retrieve the nearest line (street) for each point (sensor). Furthermore, I want to update my point-attribut "edgeid" with the id of the nearest line.

EVIS_Stat_Sensorik_.. is my point-layer and Reprojiziert is my line-layer.
For the upload string I used the "to_attr"-String because the relation between the pointlayer and the linelayer should not be based on the distance, but on the id of the nearest edge. So that the edgeid-attribut will get updated by the function and not some other attributs of the pointlayer.
For the "Column name where values specified by upload option.." I chose edgeid because the edgeid should get updated.
For the "Column name of neariest feature.." I chose id, because the attribut edgeid  of the pointlayer should be updated WITH the id of the nearest line of the line layer.

But the problem is, if I run the algorithm, the output is a NEW pointlayer which is identical with the given input pointlayer EVIS_Stat_Sensorik_.. and the edgeid-attribut of BOTH the input and output pointlayer is still empty. So I do not understand where I can see the results of the algorithm.
In the following picture you can se the attributtable of my output-linelayer named "Nearest":

I have looked through a lot threads and the duplicate questions do not help me either (How to estimate distance from points to line feature with v.distance? , Nearest distance between point layer and line layer, QGIS, Finding nearest line to point in QGIS? ). Any suggestions?


texts: Seeking testing software (like QTP, Selenium) for ArcGIS platform (Desktop and Web)?

texts:
Does anybody know any testing software/tools for ArcGIS platform (version 10.1)?

Specifically: Customized ArcGIS Desktop application & Web application.

I researched on Google but am not able to find any specific software.  Few experts told me to try manual testing so I'm a little bit confused.


texts: How to follow oneway directions using sld?

texts:
I want to show oneway directions of streets using SLD in geoserver. I follow this instruction but i still have a problem with directions of arrows. Most of all arrows get south to north and west to east direction. The arrows don't follow the direction of the lines. How to solve this issue?
<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>oneway_ir</Name>
    <UserStyle>
      <Title>A gold line style</Title>

      <FeatureTypeStyle>

       <Rule>
               <MaxScaleDenominator>2000000</MaxScaleDenominator>
      <TextSymbolizer>
        <Label>
          <ogc:Literal>&#x2192;</ogc:Literal>
        </Label>
        <Font>
          <CssParameter name="font-family">Lucida Sans</CssParameter>
          <CssParameter name="font-size">18</CssParameter>
        </Font>
        <Fill>
          <CssParameter name="fill">#000000</CssParameter>
        </Fill>
         <VendorOption name="followLine">true</VendorOption>

      </TextSymbolizer>

    </Rule>

          </FeatureTypeStyle>
        </UserStyle>
      </NamedLayer>
    </StyledLayerDescriptor>


texts:
I'm using pgsql2shp.exe from a Windows batch file and I get an error when using a WITH clause:

    "C:\Program Files\PostgreSQL\9.4\bin\pgsql2shp.exe" [ ... ] "WITH test AS (SELECT gid, the_geom FROM siveer_aep.geo_c_3_0_3) SELECT * FROM test"

The [ ... ] stands for options I consider pointless here.

I have tried a simple query that runs fine. The error output is:

Table geo_c_3_0_3) SELECT * FROM test does not exist

It looks like WITH clause syntax is misinterpreted by the exec...

I have not found anything on the Internet about this issue. Any idea?


texts: Finding and replacing strings with accented characters using ArcPy UpdateCursor?

texts:
I have a number of feature classes in which some records in the attribute tables need updating with new information.  However the strings to replace aren't found in the UpdateCursor due to the presence of accented characters.  I've run the same search on the table and just removed the accented e (é) from both the record and the search string and the string is found and replaced correctly.

# For each row, evaluate the PRESENCE field (index position
# of 0),
# and update with new value if the record reads
# "Présence n'est pas probable"

with arcpy.da.UpdateCursor(fc, inputField) as cursor:
    for row in cursor:
        if row[0] == "Présence n'est pas probable":
            row[0] = "Non probable"
            cursor.updateRow(row)
        else:
            pass

How does one search for a string with accented characters?


texts: Error Code -2147483647 when combining vector layers

texts:
I am trying to combine multiple vector layers in QGIS using the "merge vector layers" function from the dropdown lists

(vector -> Data Management Tools -> Merge Vector Layers)

All data are point data but when I run the process it returns "-2147483647 See log for more details" and when I open the log message panel it shows "key error -214783647".

Can someone please explain what this error is and how I can remedy it?

Edit
I'm using QGIS 2.18
The point data I'm using is waypoint data from 2 surveys which were uploaded to Google Earth and saved as a .kml file. I uploaded it to QGIS where each set of points shows as separate layers and would like to combine them into a single layer then save as a shapefile.
link to data:
https://drive.google.com/open?id=109K5lobpPGVHNVmG3R43nLUt1HTXHinH

I hope this clarifies some. I'm no GIS wiz.


texts: "CRS selector" dialog after running a script

texts:
I add an empty vector layer using the following script.

layer =  QgsVectorLayer('LineString', 'MyLayer' , 'memory')
QgsMapLayerRegistry.instance().addMapLayers([layer])

After I run the script, "Coordinate Reference System selector" dialog appear.  How can I get rid of this dialog?  Or, is there a script that will close the dialog?


texts: How to "lower the resolution" of a lat long csv

texts:
I have a csv with 12.5 million date/lat/long/value points, it's a point every 800 meters ( in a grid ) of the united states

what's the best way to convert the points to every 1600 meters? or 3200? etc...

Really need to reduce the file size


texts: Exporting image from ArcGIS map service where bounding box wraps?

texts:
I'm trying to export maps from an ArcGIS REST Service using WGS84 and have run into difficulty when the bounding box ends up exceeding latitude 180.

I was hoping that the map would wrap, but it doesn't. In some cases the resulting image ends up zoomed very far out and not anywhere I intended it to be.

This is an example of where the map gets clipped where I'd prefer it wrapped.

This is an example of where I've exceeded 180 and then the image is zoomed out and not on the location I intended.

Using the coordinates of a location I buffer it with a rectangle for the bounding box, which means that in some cases the bounding box will overlap -180 or +180. The resulting bounding box then doesn't get my desired result.

What is the correct way to get this to work?


texts: Landsat Surface Reflectancce product Python packages

texts: I have to work with numerous Landsat surface reflectance product that I downloaded from http://landsat.usgs.gov/CDR_LSR.php. The product comes with a QA data with information on clouds and shadows that can be used to mask the data. Has anyone here used python to apply these mask or know of any working libraries.


texts: user select fields to be deleted in a feature class

texts:
How do I make the script below accept parameter arguments correctly in ArcGIS 10? I added it to a toolbox as a script, created parameter called Fields to be Deleted: and type as Field, but when I double click it to run it, the display dows not allow for adding fields (to be deleted)? PLEASE HELP:(

# Import arcpy module
    import arcpy
    arcpy.env.workspace = "c:\\test\\Mydata.mdb"
    Drop_Field = arcpy.GetParameterAsText(0)

# Process: Delete Field (36)
    arcpy.DeleteField_management("ArcAddress", Drop_Field)


texts: Changing page settings in Layout View in ArcMap?

texts:
I need to print a map from ArcMap. I'm going to use a plotter and want to print it on a A0 format paper. But...the layout settings are totally driving me crazy!
As you can see in the first picture below there are 3 boxes I can play with:

"Printer Paper" (Which is the smallest solid one in the top right corner of the first picture)
"Printer margins" (which is the largest dashed one)
"Map Page (Page Layout)" (which is the largest solid one).

The only one which seems movable and shapeable is the Map Page. I would expect that Printer Paper and Printer Margins to match, but they don't and that seems to be the reason I can't set it right as you can see in the print preview of the second picture.

I'm sure I'm missing something. Any help?

In the picture below I unchecked "Use Printer Paper Settings" and size got inverted compared to previous situation. Still, as you can see in the print preview in the low right corner, the result is not what I want.


texts:
I want to perform a trend analysis for each pixel over it's time dimension.
Therefore I wrote a loop over each pixel doing the calculation but this is extremely slow.
I am sure there must be a faster way to perform the computation.

Also the dimensions of the file are 1440,720,480.
Is this a doable size for R? The problem is that i cannot find the statistical function implemented in IDL which I otherwise use.

gridfile = 'test.nc'
data = brick(gridfile)

require(fume)

mk = matrix(numeric(0), 720,1440)

for(i in 1:72) {
  for(j in 1:144) {
  temp = mkTrend(as.vector(data[i,j,]))
  mk[i,j] = temp$`Sen's Slope`
  }
}


texts: Secure ArcGIS connections to Access .accdb?

texts:
With ArcGIS is it possible to secure connections to an Access 2007+ format database tables?

At Connecting to a 2007 Microsoft Access database (.accdb) in ArcGIS it says "if the database has a password, enter that information in section two" so the database can be secured, but what about a table? And in the event table-level security is not needed, where are the credentials stored? Is that OLE connection just for the logged in user or does it apply to the whole machine?

If you have the requirement that Jonah Whale be allowed to read table Big_Fish_Tall_Tales while Captain Ahab can not, and ArcGIS and Access 2007+ is your toolset, what is the required and/or recommended setup?


texts: I am wondering how to or if it is possible to crop a SPOT satellite image (the image has black edges) using ArcMap 10 and if so how do you do this.


texts: Removing the Leaflet draw toolbar when added as Leaflet control

texts:
I add a draw-toolbar (I'm using the my own version from the example "Disabling a toolbar item")

drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems);
var drawControl = new L.Control.Draw({
    draw: {
        polygon: false,
        marker: false
    },
    edit: {
        featureGroup: drawnItems
    }
});
map.addControl(drawControl);

This works great but I can't remove the created toolbar with all the handlers.
How can I achieve this?

I tried:

map.removeControl(drawControl);             //dosent remove handlers
drawControl.remove();                       //dosent remove handlers
drawControl.removeToolbar();                //dosent find function
drawControl._toolbars.edit.removeToolbar(); //dosent find function
drawControl._toolbars.draw.removeToolbar(); //dosent find function
drawControl._toolbars.removeToolbar();      //dosent find function

Edit for better understanding of my problem

The documentation claims that I can add a toolbar like these:

var toolbar = L.Toolbar();
toolbar.addToolbar(map);

How ever this gives me a this.callInitHooks is not a function TypeError and it's not what i need because I want to customize the buttons depending on what the user can edit.

So I go with the drawControl version for the toolbar. But for this verion there is no function like removeToolbar().

The only thing I coud find was the first two of my "I tried" exampels. But this only remove the html container and not the class with all the handels like the removeToolbar() function. (At least thats was I think it does when I look at the code from the removeToolbar() function ).

Redraw

var drawControl = "";

function addToolbar() {
    if( drawControl !== "" ) {
        /* remove old toolbar */
        map.removeControl(drawControl);             //dosent remove handlers
        drawControl.remove();                       //dosent remove handlers
        drawControl.removeToolbar();                //dosent find function
        drawControl._toolbars.edit.removeToolbar(); //dosent find function
        drawControl._toolbars.draw.removeToolbar(); //dosent find function
        drawControl._toolbars.removeToolbar();      //dosent find function
    }

    /* draws new toolbar (with different drawnItems and e.g. marker: true)*/
    drawnItems = new L.FeatureGroup();
    map.addLayer(drawnItems);
    drawControl = new L.Control.Draw({
        draw: {
            polygon: false,
            marker: false
        },
        edit: {
            featureGroup: drawnItems
        }
    });
    map.addControl(drawControl);

    map.on(L.Draw.Event.CREATED, function( e ) {
        console.log('created');
    });
}

But because it not fully removes the old toolbar. It fires the L.Draw.Event.CREATED also from the old toolbar.

So how can I fully remove the old toolbar?


texts: Setting cell size on "Polygon to Raster" tool, in ArcGIS?

texts:
How do I properly set the cell size that I want when converting polygons to raster files?

Also on the raster file Layer Properties what exactly does the Cell Size indicate, as it doesn't seem to correlate when I measure a cell?

Here is the background to my question.

I'm converting a shapefile that represents radio signal strength to a raster file.  The shapefile is made up of polygons that each represent average signal strength inside that polygon.  Each polygon is 3 x 3 arc-seconds with a single value to represent signal strength.  In my area 3 x 3 arc-seconds equates to about a 76.8 meter x 92.6 meter rectangle.
When I convert to raster I would like to have similar resolution, so I was hoping that when I converted the shapefile each cell could be approximately the same size or smaller than the 76.8 meter x 92.6 meter rectangle polygons.

Using the "polygon to raster tool" Figure 1, it automatically sets the cell size to .0052.  What units is that in?

When I look at the new raster file, Figure 2, it clearly states that the cell size is .0052 x .0052 but I see nothing that specifies what units .0052 represents.

Measuring the individual cells however, I get measurements that don't seem to correspond to  .0052 of any units I've tried, Figure 3. Here using the pixel inspector I identified an individual cell and found it to be about 577.86 x 475.1 meter or 1,885.4 x 1,559.03 feet.

Just to be complete here are the Dataframe setting in Figure 4 showing that the units are in Feet.


texts: Getting complete attribute table in output feature class through ModelBuilder?

texts: I'm creating a tool with ModelBuilder and whenever I run it, it creates an output however the output feature class has an empty attribute table, which I need.


texts: Insert an online map source (from a website) into ArcGIS?

texts:
Is there anyway to insert an online map source into ArcGIS? There are a lot of global maps which only is accessible by their websites, and there is no save option to save a part of these maps and use in GIS software. I would like to insert an online maps into ArcGIS to work on it.
To be more clear, I am not looking for a base map. It's a map that produced and developed by an institution. I just want to import the link of this website and open this map in ArcGIS.


texts: Use QGIS to visualize OpenStreetMap data

texts: I have a osm map for a certain region. Also using osm2po I extracted the graph which includes edges and nodes into the database. Now I want to visualize the map and then visualize the shortest path. How can I visualize the map in QGIS first? I mean it doesn't directly take the osm file? How can I accomplish this?


texts: WCS band metadata description on mapserver 7,4.4 is not being outputed

texts:
I am trying to set proper metadata on a mapserver 7.4.4 instance, based on this example:

https://mapserver.org/nl_NL/ogc/wcs_format.html#grib2-output-format

And WCS paramters like:

(wcs|ows)_band_interpretation
(wcs|ows)_band_uom

That are well  described on:

https://mapserver.org/ogc/wcs_server.html#new-band-related-metadata-entries

I have implemented something like this for testing:

LAYER
        NAME "MostProbable"
        PROJECTION
            'init=epsg:4326'
        END #PROJECTION
        EXTENT -180.0 -61.06747625006139 179.9999424001524 83.99916720600001
       METADATA
            'ows_title' 'WRB most probable class'
            'wcs_label' 'MostProbable'
            'wcs_rangeset_name' 'range1'
            'wcs_rangeset_label' 'label1'
            'ows_extent' '-180.0 -61.06747625006139 179.9999424001524 83.99916720600001'
            'ows_srs' 'EPSG:4326 EPSG:152160 EPSG:3857 EPSG:54009 EPSG:54012'
            'ows_abstract' ''
            'ows_keywordlist' ''
            'wcs_formats' 'GEOTIFF_BYTE'
            "wcs_bandcount" "1"
            "wcs_band_names" "band1"

            #default values
            "wcs_band_interpretation"   "Allo Allo this is nighthawk"
            "wcs_band_uom"              "class"
            "wcs_band_definition"       "test Definition"
            "wcs_band_description"      "Come in London"
            "wcs_interval"              "0 125"
            "wcs_significant_figures"   "3"

            #specific band values
            "band1_band_interpretation" "Listen very carefully, I shall say this only once"
            "band1_band_uom"            "class2"
            "band1_band_definition"     "SpecificDefinition"
            "band1_band_description"    "This is a specific description"
            "band1_interval"            "0 255"
        END

Doing a DescribeCoverage:

https://maps.isric.org/mapserv?map=/map/wrb.map&SERVICE=WCS&VERSION=2.0.0&REQUEST=DescribeCoverage&COVERAGEID=MostProbable

I dont see any changes on the response for example uom stays the default: <swe:uom code="W.m-2.Sr-1"/>

What could be the problem for the metadata not to be displaied on the DescribeCoverage response description, maybe I am missing something ???


texts: Editing PostgreSQL/PostGIS table via join in QGIS

texts:
I have a standard PostgreSQL/PostGIS table and another table which stores a geometry. They both already have a relationship with ondelete/onupdate cascade.

I make a PostGIS connection in QGIS 2.18, I add both to the layers, I can access both of their attribute tables. I go to Project/Project Properties and set the relation.

I open the geom table and join the table. I then open the attribute table in edit mode, those columns from the 'geom table' can be edited, however I cannot edit those columns which are joined from the other table.

Where I have gone wrong?


texts: How should the raster calculator be used?

texts:
I have a set of rasters (in a geodatabase) and I'm trying to do a simple operation with them. Please see code below. I'm indicating with an arrow at the end of the code the operation I'm trying to do:

import arcpy
import os
from arcpy.sa import *

script_folder = sys.path[0]
working_folder = os.path.normpath(os.path.join(sys.path[0], '..'))

DSM_model = os.path.join(working_folder, 'Data4Optimimty\ESRI_UK_Bluesky_DSM\ESRI_UK_Bluesky_DSM.gdb', 'DSM_25cm_2015')

dsm_1 = os.path.join(working_folder, 'Optimity_Temp.gdb', 'dsm_1')
dsm_2 = os.path.join(working_folder, 'Optimity_Temp.gdb', 'dsm_2')
dsm_3 = os.path.join(working_folder, 'Optimity_Temp.gdb', 'dsm_3')
dsm_4 = os.path.join(working_folder, 'Optimity_Temp.gdb', 'dsm_4')

dsm_5 = os.path.join(working_folder, 'Optimity_Temp.gdb', 'dsm_5')
=> dsm_5 = DSM_model + dsm_3 - dsm_4

dsm_1, dsm_2, dsm_3 and dsm_4 are all existing rasters with same characteristics (extent, cell size, etc).

In ESRI help section I can find the following:

Syntax

RasterCalculator (expression, output_raster)

So I tried with the following line:

arcpy.sa.RasterCalculator (DSM_model + dsm_3 - dsm_4, dsm_5)

But that doesn't work. Error:

Traceback (most recent call last):
  File "D:\ArcGIS_Pro_Projects\Optimity_Alf\Arcpy_Scripts\Updating_Raster.py", line 42, in <module>
    arcpy.sa.RasterCalculator (DSM_model + dsm_3 - dsm_4, dsm_5)
AttributeError: module 'arcpy.sa' has no attribute 'RasterCalculator'
Failed to execute (UpdatingRaster).

Any idea on the way I should be using the RasterCalculator?


texts: QGIS weighted centroid of polygons with weights from raster layer

texts:
I am a total novice with GIS and using QGIS as that is what is available to me.

I have a vector layer of polygons (about 3500) and a raster layer of light intensity. What I'm working with looks like this:

I am trying to find the centroid of the polygons based on light intensity as a proxy for center of population in the polygon.

I know how to calculate the centroids of the polygons, but I am not sure how to include the values from the raster layer as weights.

I've read about clipping the raster layer based on a selection from the vector layer. I can clip a section from the raster and calculate the centroid based on the light intensity band, but I have no idea how to map that back to the original projection as the clipped raster only contains a GeoTransform (containing the coordinates for the top left corner of the rectangle that contains my polygon and some pixel extents), and a ProjectionRef (a string that I frankly understand little about).

I read Generate centroids for vector polygons based on raster and land-use data using QGIS, which makes sense. But looking at the interface for r.mapcalc in qgis I am at a loss of how to calculate the weighted centroid per polygon.

Logically I would assume that this kind of tool or method exists, but as I am a novice I don't know where to look and probably don't know how to ask properly either.


texts: What's the best way to implement routing through areas such as parking lots?

texts:
I'm using OSM2PO to route for walking directions and I've noticed a lot of the OSM walking paths are not ideally setup.

For example there are many paths which get divided up by parking lots, which are obviously walkable, but they are ignored and not converted to roads, so the router believes you can't walk through them.  This results in much longer walking distances than necessary.

Is there any practical way to deal with this?  I guess the same problem occurs when dealing with any "areas" that are walkable (ie. parks without paths through them).

Another issue is where sidewalks do not connect to roads.  In many situations this again leads to sub par routing results.

OSM has these issues all over the place.  It sounds like a difficult problem, but a common one, so I'm wondering if anyone has addressed this in a practical way.

Thanks for your thoughts!


texts: I want to download Geonode binaries for Windows and unfortunately its a while that all geonode.org download links for Windows doesn't work. So how can I download or find geonode-2.6.x or higher?


texts: How to eliminate the overlapped feature parts which are exposed by Union?

texts:
Here is the scenario:

I am performing Union on 2 shapefiles: A & B

A: 300 features B: 8 features.

The result of the union generated is 334. Through the helpful comments I know that it was due to some overlapping features within either shapefile.

So now my question is how to get the union results I want which will simply have 300+8= 308 features.


texts: Querying PostgreSQL with GIS to get GeoJSON for GeoJSON tile layer?

texts:
I have a PostgreSQL database with one table like the following:

CREATE TABLE post_pl
(
  gid serial NOT NULL,
  zip character varying(5),
  zip_n integer,
  city character varying(50),
  geom geometry(MultiPolygon),
  CONSTRAINT post_pl_pkey PRIMARY KEY (gid)
)

In this table there are stored polygons as geometry.

I want to get all polygons from that table within a bounding box.
I use an htaccess file to route all requests to one PHP file.
Then I use params of the url to create a bounding box.

Example:
URL: polygons.php/{z}/{x}/{y}.json

Creating bounding box(see: osm slippy tilenames):

function tileCoord2BBox($url)
{
    $parsed = explode("/",$url);

    $z = $parsed[3];
    $x = $parsed[4];
    $y = $parsed[5];

    $n = pow(2, $z);

    $lon1 = $x / $n * 360.0 - 180.0;
    $lat1 = rad2deg(atan(sinh(pi() * (1 - 2 * $y / $n))));

    $lon2 = ($x + 1) / $n * 360.0 - 180.0;
    $lat2 = rad2deg(atan(sinh(pi() * (1 - 2 * ($y + 1) / $n))));

    $bbox =  array();
    $bbox['lon1'] = $lon1;
    $bbox['lat1'] = $lat1;
    $bbox['lon2'] = $lon2;
    $bbox['lat2'] = $lat2;

    return $bbox;
}

I used those return value in ST_MakeEnvelope to limit the results.

What i want to know is:
How should a postgres statement look like to get all polygons within those
bounding box and that I can use the result as GeoJSON Tile Layer like this one: GeoJSON Tile Layer?

If possible with a JSON attribute that I can identify each polygon by "zip" on client side.

I want to have a similiar one like that one in the link above with my own geodata in Leaflet.


texts:
Is there a way to get only the relation IDs of a boundary in OSM?
E.g. relation #350303 http://www.openstreetmap.org/relation/350303 is Pärnu maakond in Estonia which is an admin_level=6 boundary.

So my question is: How do I get the relation IDs of every admin_level=X boundary in country XY?


texts:
I have a problem in the section graduate( symbology--->graduate).
I want to scaling data from 0 to 0.9 in ten step, but when i put in the cases a range from 0.0 to 0.1 (or 0.1 to 0.2 etc)  the system  transforms the range from 0 to 1.000000 (or 1.000000 to 2.000000).

I have rinstalled the program and now  I have the 3.10, but I have the same problem. I have also tested my shape file in an another computer and the approximation works. What could be the problem?


texts:
I'm adding polylines (elevation curves) to a world map. Said polylines are stored in a PostgreSQL database from gdal_contour tool. The geometry columns (LineString, 3857) has a GIST index and a CLUSTER on it. When Mapnik renders tiles, the GIST index is used.

Two bottlenecks are causing my map to be extremely slow :

Even though GIST index is used, the SELECT statement takes several seconds (usually around 30 seconds) to fetch results.
When the SELECT statement is done, renderd takes even more time to render the tiles (around 90 seconds).

My first guess is, the contour lines are quite long, and for most of these, a big part of the geometry is probably outside of the tile. I'm guessing Mapnik works with the entire shape, even if only a small part should be drawn on the tile. If I'm right, how should I split (or import) my contour shapes correctly, to work with smaller parts instead? (I need to draw text along the polyline with text-spacing at 500, so lines should not be shorter than this).

My second guess is, since the table is quite big (currently 45 341 726 rows), the GIST index is not enough (but could be sufficient if the polylines are splitted maybe?). I'm not sure how to handle this, but I'm wondering if there's a way to split the index in several parts, and index the splitted index to speed things up? Or is this a case when GeoHash would be helpful?

The polyline layer is defined as follows :

- id: contour100
    geometry: linestring
    Datasource:
      type: "postgis"
      dbname: "gis"
      key_field: ""
      geometry_field: "wkb_geometry"
      extent: "-20037508,-20037508,20037508,20037508"
      table: (SELECT
            wkb_geometry, height
            FROM contour100
            WHERE wkb_geometry && !bbox!
          ) AS c100
    properties:
      minzoom: 12

The server has 128 GB of RAM, and Xeon E5-1660v3 processor. The database is stored on NVMe drives, so the problem has to come from a wrong configuration on my side.


texts: PostGIS : Update a POINT (geometry) with a POINT from another table with the same SRID

texts:
I try to update a geometry point with data from another table (also a geometry point).
Both are SRID 4326.
The query is executing without errors, everything is updated but the geometry is just not updated.
Do I have to transform first to lat lng and then back again?
Is it not possible to use a geometry type?
I don't get it at the moment.
Example:
UPDATE location
    SET
    location_name = import_ogr.location_name,
    location_category_id = import_ogr.location_category_id::INTEGER,
    position_geo = import_ogr.wkb_geometry
    FROM importexport.import_ogr
    WHERE location.location_id = 23017;


texts: Masking using two sliders on Google Earth Engine

texts:
I am trying to apply a mask based on two conditions (elevation and MNDWI), however, I don't find a way to integrate the two conditions using sliders, which would update the map when any of the sliders are used.
I have created the sliders and variables however they work independently, applying only one of the two thresholds when a slider is used.
At the end of the code, the variable "desired" describes the desired output.

    var start = '2014-09-15';
    var end = '2014-12-30';
    var year = '2014_wet';

    var ndwicollection = ee.ImageCollection('LANDSAT/LC8_L1T_8DAY_NDWI')
                         .filterDate(start, end).first();
    print ('ndwicollection', ndwicollection);

    var srtm = ee.Image('USGS/SRTMGL1_003');

    var mndwiParams = {min: -1, max: 1, palette: ['green', 'yellow', 'blue']};

    // threshhold

    // Create a label and slider for MNDWI.
    var label = ui.Label('MNDWI Threshold');

    var slider = ui.Slider({
      min:-1,
      max: 1,
      step: 0.01,
      style: {stretch: 'horizontal', width:'500px' },
      onChange: (function(value) {
      Map.layers().reset();
      var ndwi_collection_mask =  ndwicollection.updateMask
      (ndwicollection.gte(value));
      Map.addLayer(ndwi_collection_mask, {palette: 'blue'}, 'MNDWI threshold '+String(value));
      return ndwi_collection_mask;})
      });
      slider.setValue(0.0);  // Set a default value.

    // Create a label and slider for Elevation.
    var label1 = ui.Label('Elevation');
    var slider1 = ui.Slider({
      min:0,
      max: 3000,
      step: 1,
      style: {stretch: 'horizontal', width:'500px' },
      onChange: (function(value1) {
      Map.layers().reset();
      var srtmmask =  srtm.updateMask(srtm.lte(value1));
      Map.addLayer(srtmmask, {palette: 'blue'}, 'MNDWI threshold '+String(value1));
      return srtmmask;})
      });
      slider1.setValue(0);  // Set a default value.

    // Create a panel that contains both the slider and the label.
    var panel = ui.Panel({
      widgets: [label, slider, label1,slider1],
      layout: ui.Panel.Layout.flow('vertical'),
      style: {position: 'bottom-right',width: '520px'
          }
    });

    // Add the panel to the map.
    Map.add(panel);

//// to replace with sliders
    var desired = srtm.updateMask(srtm.lte(100)).and(ndwicollection.updateMask
    (ndwicollection.gte(0)));
    Map.addLayer(desired, {palette: 'blue'});


texts: ArcMap cannot access WCS version 2.0.1 on MapServer?

texts:
I have a MapServer 7.0.7 map containing only a layer with a GeoTiff DEM.

It works perfectly with QGIS 3.0.3 but I cannot use it with ArcMap 10.5.

This is my map file:

MAP
  OUTPUTFORMAT
      NAME "GTiff"
      DRIVER GDAL/GTiff
      MIMETYPE "image/tiff"
      IMAGEMODE FLOAT32
      EXTENSION "tif"
    TRANSPARENT ON
  END

  NAME "Modèle numérique de terrain"
    EXTENT  -6.0001389 39.9998611 10.0001389 52.0001389
    MAXSIZE 4096
    PROJECTION
        "init=epsg:2154"
  END

  WEB
    METADATA
            "wcs_label"           "WCS Server"
            "wcs_onlineresource"  "http://localhost/cgi-bin/mapserv?map=/tmp/testNaN.map&"
            "wcs_enable_request" "*"
            "wcs_fees"            "none"
            "wcs_accessconstraints"    "none"
            "wcs_keywordlist"          "wcs,dem"
            "wcs_metadatalink_type"    "TC211"
            "wcs_metadatalink_format"  "text/plain"
            "wcs_metadatalink_href"    "http://localhost"
            "wcs_address"              ""
            "wcs_city"                 "Paris"
            "wcs_stateorprovince"      "IDF"
            "wcs_postcode"             "75000"
            "wcs_country"              "France"
            "wcs_contactelectronicmailaddress" "admin@gnone.world"
            "wcs_contactperson"            "me"
            "wcs_contactorganization"      "unemployed"
            "wcs_contactposition"          "manager"
            "wcs_contactvoicetelephone"    ""
            "wcs_contactfacimiletelephone" ""
      END
    END

  LAYER
        NAME TESTNAN
          DATA "/tmp/tmp.tif"
        STATUS ON
        TYPE RASTER
        PROCESSING "SCALE=AUTO"
        PROJECTION
           "init=epsg:2154"
        END
    METADATA
      "Altitude_band_description" "Altitude en mètres"
      "wcs_band_names" "Altitude"
      "wcs_bandcount" "1"
      "wcs_description" "Modèle numérique de terrain"
      "wcs_enable_request" "*"
      "wcs_label" "IGN BDAlti V2 (rés: 75m)"
      "wcs_name" "BDALTI75"
      "wcs_rangeset_label" "IGN BDAlti V2 (rés: 75m)"
      "wcs_rangeset_name" "BDALTI75"
      "wcs_rangeset_nullvalue" "-99999"
      "wcs_srs" "EPSG:2154"
    END
  END
END

While I add my WCS server under ArcMap 10.5 without forcing the version. ArcMap makes two HTTP GET request

:/1 - - [14/Jun/2018:16:56:14 +0200] "GET /cgi-bin/mapserv?map=/tmp/testNaN.map&SERVICE=WCS&SERVICE=WCS&REQUEST=GetCapabilities HTTP/1.1" 200 6816
::1 - - [14/Jun/2018:16:56:14 +0200] "GET /cgi-bin/mapserv?map=/tmp/testNaN.map&SERVICE=WCS&SERVICE=WCS&VERSION=2.0.1&REQUEST=DescribeCoverage HTTP/1.1" 400 536

As you can see the second request does not work and Mapserver answer a HTTP 400 error.

While entering this request in Chrome I can have the full error:

msWCSDescribeCoverage20(): WCS server error. Missing COVERAGEID parameter.

the only workaround I found is to force the version to 1.0.0

How do I make MapServer understand ArcMap WCS 2.0.1 protocol?


texts: Draw arrow in direction of line on midpoint in Geoserver SLD

texts:
Similar to this question: How do I draw the arrow on the line using SLD of Geoserver?

The answer there works, but I need to place the arrow on the midpoint of the line. Alternatively, it could be placed along the line with a certain offset from the start or end point.

The reason is that I am using this for lines that have a symbolized point feature (for manholes in a sewer network) at each end. These points would interfere with arrow at beginning or end of the line.


texts: Applying osm2pgsql .diff patch in windows 7

texts:
When attempting to run osm2pgsql in command prompt, I ran into the password prompt as described here. This requires installing a patch. I'm unfamiliar on how to install a .diff patch on Windows 7.
Do I need to install GnuWin as shown here?
Thanks for all the help!

EDIT: This is the same question as this one on Stack Overflow.


texts: I need free terrain data for the whole world which I can import into TileMill. Is there some source where I can download them?


texts: Limiting list values based on previous field in QGIS Forms?

texts:
Looking for assistance to set up appropriate configuration to enable drop-down lists values to be filtered based on a previous value.
Basically, I have a 3-tiered categorisation structure to one of our GIS layers - Class > Type > Sub Type.
Once Class has been selected I would like the Type list to be filtered to only values relevant for that class and so on for Sub-Type.

I have turned this into Code Value list attribute tables and loaded inot QGIS as attribute tables -

I have used these tables as Value Relation lists in the widget properties

The drop down lists work fine without a filter. But as soon as I try adding a filter with an expression, there is nothing in my drop down...
If I simply type "CLASS_ID" = 1 ... I do limit the list to only those for Class = 1. The expression on it's own (attribute($currentfeature, 'LU_CLASS') correctly populates the integer values if I use it to populate a virtual field in the expression builder...
Can anyone provide any guidance as to why this doesn't work in the widget filter?
I have also (briefly) tried to set up a relation in the project properties... but this feels like I am going about it the wrong way... the records aren't a child or related to this other table - I simply want to use them to populate drop-down lists.
Am I just missing a simple step?


texts:
I have a question with regard to spatial aggregation in R.
What I am trying to do is aggregate a point dataset to a grid.
I am unsure however how to do this as I have little experience with this sort of stuff.
I was hoping anyone of you might have some useful guidance/a possible solution.

My vantage point is a dataset containing georeferenced data on conflict events in Africa (see www.acleddata.com).
The points are georeferenced with latitude/longitude coordinates and contain data on event type and time.
What I want to do is aggregate these points to a 1x1 degree grid.

Thus a grid-cell should contain the information of the data points if a event happened to occur within that grid-cell. The eventual product of this should be a data frame or something that I can export to a csv-file as the data is intended to be used in a panel data-set for statistical analysis.

So far I loaded and plotted the data and the shapefile using the code below. I believe that I should use the over function from the sp package to
aggregate but I do not know how. Hope one of you can help.

The code I used so far can be found here with the corresponding visual result over there.

Suggestion for doing this in QGIS are welcome as well.


texts:
I have difficulty copying the FID column known as OBJECTID to another column.

arcpy.CalculateField_management("file.shp", "index_field", "OBJECTID")

if this run it gives a window that says it can't copy the content. 'The calculated value is invalid for the row with ObjectID = 1'. What is the problem?

it gives this error for any column not just FID.


texts: Why does my reprojected globe visualisation not draw?

texts: I have followed all of the instructions in QGIS Map Design, Graser and Peterson, Locate Press 2016, and cannot get the visual of the reprojected shapefiles clipped to Rome coordinates to draw as in Fig 12.4 p86.  Is there another step beyond applying the CRS as describe on p85 that I need to do? I have used 2.8, 2.14, 2.16 without success.


texts: How do I reclassify raster pixels that are not intersected by a line?

texts:
I have a large raster and a line.

I want to preserve the values of the pixels that are below the line, and reclassify all other pixels to NO VALUE or 9999, so that only the pixels where the line overlapped are the same.


texts: Does exist a Python library to generate polygons avoiding antimeridian issues?

texts: If I need to generate a polygon that crosses the anti-meridian, does it exists a Python library that treat such cases (partitioning such polygon and generate a multipart polygon)?


texts:
I am trying to compile QOCI in QGIS master, I followed the instructions, directly from QT, but I still get the same error: "fatal error: QtSql / private / qsqlcachedresult_p.h".
I am working on a 32 bit machine, linux Xenial and with the instant client of oracle version 12.2,

How can I fix the error?


texts: Suggestions for an efficient way to store raster files for lab use?

texts:
I'm trying to create a repository for our raster files from various sources and Im exploring the possibilities on how to store them and deliver them in our local intranet. I have come across rasdaman and seems like the tool for that kind of job.

Any suggestions on how to create a system like that would be really helpful

I'd prefer open source solutions over proprietary ones


texts: Counting self intersection in same shapefile with ArcMap

texts:
Is it possible to count the number of intersections in the same linear shapefile with ArcMap 10.2?

My license level is advanced.


texts: QGIS custom script to run one process or another depending on input

texts:
[QGIS 3.10] I have two processes I have created using the graphic modeler.
I'm trying to create a script that would run one process or another depending on the inputs. But I'm getting an error in line 30 regarding the process.

1) Is this because it is not recognising the models I've created? any idea how can I solve this?

2) Also, any idea how can I restrict the Vector type of the Inputs?

  NameError: name 'processing' is not defined

EDIT: As @Kadir suggested, I added import processing at the beginning and it seems to be working. I'm still not sure how to restrict the Vector_Layer to points or lines when using @alg decorator, in case someone can help. Else, I'm closing this question.

from qgis import processing
from qgis.processing import alg
from qgis.core import QgsFeature, QgsWkbTypes, QgsField, QgsFields, QgsPointXY, QgsGeometry, \
QgsCoordinateReferenceSystem, QgsProcessingMultiStepFeedback, QgsVectorFileWriter, QgsVectorLayer, \
QgsProcessing,QgsProcessingAlgorithm,QgsProcessingParameterFeatureSink
from PyQt5.QtCore import QVariant

@alg(name="Test if", label=alg.tr("Test if"), group="scripts", group_label=alg.tr("Example Scripts"))
@alg.input(type=alg.VECTOR_LAYER , name="LINES", label="Interpretation Lines", optional=False)
@alg.input(type=alg.VECTOR_LAYER , name="POINTS", label="Interpretation Points", optional=False)
@alg.input(type=alg.VECTOR_LAYER , name="KP", label="KP", optional=False)
@alg.input(type=alg.VECTOR_LAYER , name="DOB", label="DOB", optional=True)
@alg.input(type=alg.VECTOR_LAYER_DEST, name="OUTPUT", label="Output layer")

def testalg(instance, parameters, context, feedback, inputs):
    """
    Description goes here. (Don't delete this! Removing this comment will cause errors.)
    """
    feedback = QgsProcessingMultiStepFeedback(1, feedback)
    results = {}
    outputs = {}

    if parameters['DOB'] is None:
       # PINS_no_DOB
        alg_params = {
            'interpretationlines': parameters['LINES'],
            'interpretationpoints': parameters['POINTS'],
            'kppoints': parameters['KP'],
            'qgis:refactorfields_1:test': parameters['OUTPUT']
        }
        outputs['OUTPUT'] = processing.run('model:FGBML_PINS_no_DOB', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['OUTPUT'] = outputs['OUTPUT']['qgis:refactorfields_1:test']
    else:
        # PINS_&_DOB
        alg_params = {
            'dob': parameters['DOB'],
            'interpretationlines': parameters['LINES'],
            'interpretationpoints': parameters['POINTS'],
            'kppoints': parameters['KP'],
            'qgis:refactorfields_1:test': parameters['OUTPUT']
        }
        outputs['OUTPUT'] = processing.run('model:FGBML_PINS_&_DOB', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['OUTPUT'] = outputs['OUTPUT']['qgis:refactorfields_1:test']

    return results


texts: Line intersection resulting in offset point


texts:
I have a co-worker who manages a large MS Access database that contains X & Y fields for various locations he tracks.  We have previously exported the data to CSV and used the QGIS Delimited Text Plugin to plot the locations.

We would now like to be plot this data directly in QGIS by connecting to the MS Access db via an ODBC connection (as opposed to export to CSV).

Can QGIS consume the data from the ODBC connection to plot the X-Y locations?
Can it be "live linked" to the ODBC connection so that when a project is re-opened, the layer "replots" all old and newly added locations dynamically (as opposed to exporting to shapefile, which would create a snapshot in time)?
Can you create tabular ODBC connections that can be joined to spatial tables within QGIS?

You can do all of this in ArcGIS ("Add XY Data" tool), so it would be great if this capability existed in QGIS. I just cannot find anywhere that mentions its existence.

I am not planning to move to a spatial database.  As a GIS person myself, yes I know that it makes sense to move to a spatial database, and yes I have used PostGIS to manage this kind of data.  But these are not options right now.  The user is not a GIS person nor wants to be, and he manages his own database in MS Access and that will not change.  I am and am pulling him kicking and screaming into the GIS world and want to make it easy as possible.  He can live with exporting to CSV, I was just hoping there is a better solution USING ODBC FOR MS ACCESS and displaying the XY data in QGIS.


texts:
Does anyone know where I could find an up to date data set of amenities for the UK?

I'm working on a study which needs to understand the distributional impact of a potential development proposal, in this instance in the Greater Manchester area.  As part of this I need to obtain local amenity information including but not exclusive to the following:

Doctors/GPs
Schools
Nurseries
Playgrounds
Parks/Open Spaces
Hospitals
Care Homes/ Day Centres
Community Centres

From previous GIS projects I was aware of the Department for Transport core accessibility indicator data set but the latest set dates back to 2010 which is quite dated.  Also I've looked at OSM data, looking at building type and extracting the types I'm interested in but I think the information is very limited as apparently there isn't a single doctors within the Manchester Salford Inner Relief Route which I know isn't correct.


texts:
It is said in manual, that

  You can add pictures to your layout. The insert picture command can be
  found on the Insert tab in the Graphics group.

Unfortunately, I don't have any Graphics group in insert tab:

How to accomplish?


texts: Good practices for storing connection info while deploying geoprocessing tool on Server

texts:
I have a geoprocessing tool written in python.

It will be deployed on ArcGIS for Server 10.1 as a GeoProcessing Service and it's processing workflow includes some manipulations with data in SDE geodatabase.

Right now I've decided to save .sde file of connection (with login and password saved in it) and deploy it together with my script to ArcGIS for Server.

Another similar option is to deploy a .lyr file. In any of these cases deployment on another server requires additionally editing these files in ArcGIS for Desktop.

Can you advice any other methods for sharing data access to SDE GDB without exposing sensitive information (password)?


texts: Offseting lines with same symbology in ArcGIS Desktop?

texts: I'm working in ArcMap 10 and have been tossed an electrical system map drawn in CAD to recreate in ArcMap. In several sections, different line numbers share the same structures for long stretches. Many times these have the same voltage (which is what they are symbolized on). Is there a way for me to continue symbolizing these lines by voltage while showing 3 separate lines that are not stacked for one section, 2 separate for another, etc?


texts: How to (quickly) count the number of rows in ArcGIS attribute table?

texts:
How can I quickly see how many rows of data are in my attribute table in a large file without having to scroll right down to the bottom?

Is there a way of removing the * so that I can quickly see total?


texts: proj.4 — problems to understanding the projection result

texts:
I am using pro4js for a webproject to project data obtained from the Openstreetmap database. Basically I have this code (javascript/ES 6):

const
  //some coordinates
  we = {
    lon: 50.9802813,
    lat: 11.3290855
  },

  //the projectors
  { forward,
    reverse } = proj4('EPSG:3857')
;

let
  [  lon1,
     lat1  ] = forward([we.lon, we.lat]),

  [  lon2,
     lat2  ] = forward([we.lon + 1, we.lat + 1])
;

lon2 - lon1 = 111319.49079327378
lat2 - lat1 = 113736.3960447507

Why are those values this huge? What unit is that? Since these are projected coordinates they should represent coordinates on a flat 2d plane, don't they?


texts: Unable to execute Processing nativecentroid in python standalone scripts using PyQGIS

texts:
I am using QGIS 3.2 on Windows7. I had successfully set all enviroment variables in OSGeo4W shell. I have run small scripts without any issue.
I have to run

processing.run('native:centroids', {'INPUT':'C:/Test/demand/valid.shp','ALL_PARTS':True,'OUTPUT':'C:/Test/demand/centroid.shp'})

The above line works perfectly in console within QGIS. But when I try to run in standalone scripts it is giving following error:

Although, in my script I have added below lines and I was able to execute processing.run('qgis:voronoipolygons', parameters)

import processing
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

I don't understand where I am doing wrong.


texts: QGIS Measure Distance Using Field Calculator

texts:
I have a point shapefile with 80K+ records.  Shapefile is projected in wgs84.  The shapefile contains fields with the latitude and longitude of the point.  I need a new field that contains a distance from a known point (Point B) in each record.  Point B is the same for every record in the shapefile.  There is a distance function in the field calculator (field calculator -> geometry -> distance).  But I can't make sense of the syntax.

Example:

distance( geom_from_wkt( 'POINT(4 4)' ), geom_from_wkt( 'POINT(4 8)' ) ) → 4

What I want to do is measure the distance from the geometry of the current point to the geometry of a point I supply (for example point at "25.906543, -80.546800").  How do I refer to the current record using that tool and how to a create the new geometry?  Is is something like:

distance( geom_from_wkt( 'POINT(25.906543, -80.546800)' ),
($geometry ) )


texts: What is geometry(4326), and can you do Cartesian calculations with unprojected ellipsoidal coordinates?

texts:
I know Geography(4326) is an unbounded, and unprojected, and covers the entire world with ellipsoidal coordinates. I get that the math on this system is done on a spheroid. I also get that with a geometry the math is done on a Cartesian grid. What I really would like to understand is how SRS 4326 can be subject to both ellipsoidal calculations, and Cartesian calculations without reprojection  (like ST_Transform)?

SELECT ST_Distance( ST_MakePoint(0,0)::geography, ST_MakePoint(1,1)::geography );
   st_distance
-----------------
 156899.56829134
(1 row)

-- same result as ST_MakePoint(x,y)::geography::geometry
-- same result as ST_SetSRID( ST_MakePoint(x,y)::geometry, 4326 )
-- same result as ST_SetSRID( ST_MakePoint(x,y), 4326 )
SELECT ST_Distance( ST_MakePoint(0,0)::geometry, ST_MakePoint(1,1)::geometry );
   st_distance
-----------------
 1.4142135623731
(1 row)

By extension, MySQL does not support nearly the amount of Geographic operations as PostGIS. Worse, it determines the type of calculation (geometric vs geographic) based not the type (as with PostGIS's geometry vs geography) but based merely as the SRS ID marker. For any operations that do not support a geodetic reference system, you have to use a Cartesian reference system (which is obviously major suck because MySQL can't reproject). But why if it has these two shortcomings (fewer geodetic calculations), and an inability to reproject can it not do a geometric calculation (Cartesian reference system) with SRS ID 4326? Lastly, is there any SRS ID that is not geodetic that behaves like PostgreSQL geometry type in SRS ID 4326, so I can use it with MySQL?


texts: Counting points in polygon using weight and class simultaneously

texts:
I am using QGIS 3.1. I am wondering if there is a way to count points in polygons but by using both the weight and class.
For example, I have distributions of a number of species, each with a unique 'score'. I want to count each occurrence of a unique species (class - species), then weight each occurrence by its score (weight - score) in a grid.
Doing each of these is very simple using the count points in polygon feature but I can't do both at once.


texts:
I had drafted a polyline with QGIS, I want to continue that line with snapping from it's end, but snapping is not working in starting, after starting I am getting snapping.

Using windowsXp and Q gis 1.8.0

My issue is I had drafted a line and finished entering all attributes then i wan't to start a new line from the end point of first one there i didn't getting snapping.If i had started drafting with out snapping, second pic i can get snapping where i first wish to snap.


texts: Openlayers getExtent() returns wrong values when the dataset crosses the dateline

texts:
I have a geojson with some data which I loaded into Openlayers.
When I try to get the getExtent of the data source I get awkward values not corresponding to the extent of the data.
The following image shows the dataset in blue and the extent points, calculated with the gjsonLayer.getSource().getExtent() in red.

When I click in the polygons in the Chatham islands (on the right side of the image) I get the coordinates [183.49365234375, -43.905808356157394], out of the range -180, 180. Probably that's the reason why the method getExtent() returns messed up values.

Just to show that the problem is the data crossing the dateline, below is an image showing the dataset without those polygons on the right (in the Chatham islands) that pass the dateline. The red points are the result of the getExtent() method and they are correct.

Any idea how to make the getExtent() method return correct values when using a dataset in which some data passes the dateline?


texts: Exporting time series of multiple polygons to .csv. format table

texts:
I wanted to export a time series of Rainfall (CHIRPS) on *.CSV format without using the ui.Chart.image.seriesByRegion command. I cannot build the chart because there is too much data (daily data from 2012-2017) and I get an error.

Similar to what have been proposed here.

I'm using the following code, from which I get in rows the date, and by column the value corresponding to my polygon. However in the last column (.geo) it is the coordinates of each polygon.

Because my polygons are complex (not as in this example below) the document is becoming very heavy (around 1 Gb).

How do I remove this last column before exporting the data?

var regions = ee.FeatureCollection([
  ee.Feature(    // San Francisco.
    ee.Geometry.Rectangle(-122.45, 37.74, -122.4, 37.8), {label: 'City'}),
  ee.Feature(  // Tahoe National Forest.
    ee.Geometry.Rectangle(-121, 39.4, -120.8, 39.8), {label: 'Forest'}),
  ee.Feature(  // Black Rock Desert.
    ee.Geometry.Rectangle(-119.15, 40.8, -119, 41), {label: 'Desert'})
]);

var temps2013 = ee.ImageCollection('LANDSAT/LC8_L1T_32DAY_TOA')
    .filterDate('2012-12-25', '2013-12-25')
    .select('B11');
var triplets = temps2013.map(function(image) {
  return image.select('B11').reduceRegions({
    collection: regions.select(['label']),
    reducer: ee.Reducer.mean(),
    scale: 30
  }).filter(ee.Filter.neq('mean', null))
    .map(function(f) {
      return f.set('imageId', image.id());
    });
}).flatten();
print(triplets.first());
var format = function(table, rowId, colId) {
  var rows = table.distinct(rowId);

  var joined = ee.Join.saveAll('matches').apply({
    primary: rows,
    secondary: table,
    condition: ee.Filter.equals({
      leftField: rowId,
      rightField: rowId
    })
  });

  return joined.map(function(row) {
      var values = ee.List(row.get('matches'))

        .map(function(feature) {
          feature = ee.Feature(feature);
          return [feature.get(colId), feature.get('mean')];
        });
      return row.select([rowId]).set(ee.Dictionary(values.flatten()));
    });
};

var table1 = format(triplets, 'imageId', 'label');

var desc1 = 'table_demo_';
Export.table.toDrive({
  collection: table1,
  description: desc1,
  fileNamePrefix: desc1,
  fileFormat: 'CSV'
});


texts: Adding extra parameter in QGIS processing in grass to close files after processing

texts:
The code below works but as this involves averaging thousands of files, grass shows an error message:

Too many open files

The code:

ymax=31.4648437503
ymin=17.4068080185
xmax=-80.4330407378
xmin=-98.0539218184
raster_list = create_qgis_raster(files)
general.runalg('grass7:r.series', raster_list, False, 0,
   '-10000000000,10000000000',
   '%f,%f,%f,%f' % (xmin, xmax, ymin, ymax),
   0, output_path
)

The error happens because grass don't close these raster layers. The documentation of grass r.series algorithm has a flag called -z that means Do not keep files open. However, there is no place to put this option in runlang method of processing.

When I add the parameter -z, I get the error

Error: Wrong number of parameters
ALGORITHM: r.series - Makes each output cell value a function of the
values assigned to the corresponding cells in the input raster layers.
   input <ParameterMultipleInput>
   -n <ParameterBoolean>
   method <ParameterSelection>
   range <ParameterString>
   GRASS_REGION_PARAMETER <ParameterExtent>
   GRASS_REGION_CELLSIZE_PARAMETER <ParameterNumber>
   output <OutputRaster>

QGIS version 2.18

Is there any way of adding those additional flags of r.series in processing runlang method?

Edit:
Joseph's answer helped in adding the extra parameter but the added parameter has not closed open files after processing.


texts:
What does this expression mean? I see this in many python scripts.

try:
    ...
except Exception as e:
    import traceback
    map(arcpy.AddError, traceback.format_exc().split("\n"))
    arcpy.AddError(str(e))

traceback.format_exc().split("\n") returns ['None','']


texts: How to check out arcgis 10.2 desktop license using arcpy?

texts:
Referring to ArcGIS help  it seems possible to check out advanced desktop license. But typing
arcpy.CheckOutExtension('ArcInfo')

gives

Unrecognized extension code ArcInfo


texts: if column name = this, then search for a row >0 and change symbol, in QGIS

texts:
Excuse me as I am new to GIS. I am using QGIS to do two things. Firstly, from a csv master file, I have plotted points on a map based on their geo-reference and graduated their colour based on the level of species diversity. Cool.

The next thing I wish to do is change the symbology based on the presence of certain species. The data is formatted as follows, each row is a 10m x 10m plot in the field and each column is a species name (500+ columns). Where a species is present in a plot, the value is > 0.

I would like to use a list of species names located in another file to search through each column name (field label?), and where the name matches, search each row/plot for a value > 0. If the column contains a value > 0, the symbol of the point corresponding to that row would change. I.e. If a species is found in a plot, the symbol changes to indicate its presence.

I imagine this is quite a straight-forward task, however, I cannot find any information or guides for how one might achieve this. This is likely due to my poor understanding of GIS terminology. While I am a competent coder in programs such as R, this is my first forray into GIS.

I have acheived what I wanted by first processing and formatting the data in R using this code:

setwd("~/Uni/Data/")
CMRA_feb2014_uniform <- read.csv("~/Uni/Data/CMRA_feb2014_uniform.csv", stringsAsFactors=FALSE)
OrchidNVSNames <- read.csv("~/Uni/Data/OrchidNVSNames.csv", stringsAsFactors=FALSE)
CMRA_feb2014_Orchids = data.frame(NULL,NULL,NULL)
SpeciesCount = 0
i=1
j=1
k=1
x = 1
for(i in 1:length(CMRA_feb2014_uniform)){
    for(j in 1:length(OrchidNVSNames$NVSCode)){
        if(identical(names(CMRA_feb2014_uniform[i]),OrchidNVSNames[j,1]))  {
            for(k in 1:nrow(CMRA_feb2014_uniform)){
                if(CMRA_feb2014_uniform[k,i] > 0){
                    CMRA_feb2014_Orchids[x,1] <- as.integer(CMRA_feb2014_uniform[k,"Plot"]) #Plot number
                    CMRA_feb2014_Orchids[x,2] <- OrchidNVSNames[j,"NVSCode"]
                    CMRA_feb2014_Orchids[x,3] <- OrchidNVSNames[j,"SpeciesName"]
                    CMRA_feb2014_Orchids[x,4] <- CMRA_feb2014_uniform[k,i]
                    for(l in 1:length(CMRA_feb2014_uniform)){
                        CMRA_feb2014_Orchids[x,4+l] <- CMRA_feb2014_uniform[k,l]
                     }
                x = x+1
                }
            }
        }
    }
}
colnames(CMRA_feb2014_Orchids)<- c("OrchidPlot","OrchidNVSCode","OrchidName","Orchid%C",c(names(CMRA_feb2014_uniform)))
write.csv(CMRA_feb2014_Orchids, file = "CMRA_feb2014_OrchidData.csv", row.names = FALSE)

I then imported the final cvs as a layer, and styled the points using point displacement.

While I have achieved what I wanted, it would be great to know if there is a way to do this within QGIS, or be pointed in a direction where I may gain the background to carry out such tasks


texts:
I have been facing the issue of transform coordinates from a custom spatial reference to a well known spatial reference (let's say WGS84). To do that using ArcObjects is quite straightforward, but I am trying to do it without that library.

I must use some library designed for .NET so making a little research I found DotSpatial.Projections and ProjNet.

Both of these open source alternatives lack extensive documentation, so some effort is required to understand the process of implementation, although is not all that difficult to.

My question is: In your experience, Which one is better? taking into account ease of use, performance, maintainability (of the library), and support.

Any other suggestion is welcome.


texts: Unable to attach a projection before creating a GDAL image

texts:
I'm unable to add a projection to a newly-created raster with GDAL-python, using the code below:

arr = np.random.rand(600, 400)
filename = 'path_to_file_'
driver = gdal.GetDriverByName('GTiff')
x_size = arr.shape[1]
y_size = arr.shape[0]

dataset = driver.Create(filename, x_size, y_size,
                        eType=gdal.GDT_Float32)
dataset.SetProjection(projection)
dataset.GetRasterBand(1).WriteArray(arr)

It's expected in the code above that the projection isn't attached to the dataset, because the dataset has already been written when the projection is added to it.
Is there another way to create an image in a way that the projection is attached to the saved image?

Note that the projection (WGS84) that I'm trying to attach to the image to create was extracted from an image used as input to create the newly-created image (i.e. I'm just trying to keep the projection of the source image in the image to save). Here is the string returned by dataset.GetProjection() (i.e. the content of the variable projection):

  'GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS
  84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433],AXIS["Geodetic
  longitude",EAST],AXIS["Geodetic
  latitude",NORTH],AUTHORITY["EPSG","4326"]]'


texts: Are there any tools that can identify ravines/steep elevation drops from contour line data?

texts:
I am trying to create polygons around rivers that extend all the way up to the top of the surrounding slopes.

I am currently able to do this using stream data and contour data, but I would have to digitize it by hand. Because this would be a very time consuming process, I was wondering if there was an automated process of identifying ridges within ArcGIS or QGIS. I have attached two figures to show what I mean.

In Figure 1, you can see the contour lines surrounding the river. From this, it is very clear from the contours that there is a steep slope on either side of the river.

In Figure 2, I have filled in the area with red showing the desired polygon I would like to make.

Does anyone know if the software can identify these areas of steep elevation change using contour data and create a polygon of the area between them? Or, would my only choice be to digitize?


texts: Zonal Statistics as table produces wrong result (only one row of result)

texts:
I have temperature as a raster file and different ZIP Code Tabulation Areas (ZCTAs) as polygons. I'd like to calculate the mean temperature per zip code.

I need to know the temperature value per ZIP code region and I am trying to use the zonal statistics as table but my results are showing only one row out of 1769 ZIP codes. What could be the reason why this is happening? I also tried to do zonal statistics and wasn't successful either.

I thought that maybe the problem was that the raster (TIFF) is not only for California and my shapefile is of CA. I don't know if that's relevant.


texts:
I am trying to dissolve grid spaces by category using GeoPandas' dissolve function, but it doesn't finish after hours of running and I want to make sure I'm not doing something wrong and/or there isn't some way to speed it up.

Here's what I've got so far:

outputData = gp.GeoDataFrame(outputData[['category','geometry']])
outputData['geometry'] = outputData['geometry'].apply(wkt.loads)
outputDataFrame = outputData.dissolve(by='category')

Here's the breakdown: I've got a Dataframe called outputData with many variables/columns, so I limit it to the two I need and convert it into a geoDataFrame.   In order to use the geometry column, I need to apply wtk.loads else I get an error about it being a str type.  There are 12 values in the category column (0-11) and around 200,000 rows in the DataFrame.

I've run dissolve before (e.g. to get region boundaries from counties) but never on so many geometries.  I'm thinking that because in this case the geometries are all simple and non-overlapping (a grid) there may be a better way to combine them.  But if dissolve is just groupby + unary_union, then I don't know what could be better than that.

Is this a known unsolvable problem with large datasets or are there any tricks to speed it up significantly (e.g., by setting an option or using a different library/function)?


texts: Creating flow map with matrix table using QGIS?

texts:
I have a map http://prntscr.com/95zoyh and also data in a table (matrix table ) and I would like to have a map with different weight arrows depending of the data.

Is that possible to do using QGIS?


texts:
I have an ArcGIS Online subscription, however I continue to get an error similar to the one below when attempting to geocode a dbf of street addresses in ArcGIS for Desktop. I tried both 10.3 & 10.2.2.

  ArcGIS Online subscription is required for geocoding the table of
  addresses

I was able to geocode 11k addresses at one point yesterday, however I have only gotten this error in subsequent attempts - any ideas?


texts: SQL and CARTO - Dot Density Visualization

texts:
I need help with running a SQL on CARTO.

I want to generate a dot density similar to the link here.

The second portion of the SQL code deletes the cartodb column and nothing changes on the map.

SELECT
  ST_TRANSFORM(cynthia.by_state_population(the_geom,DP0080005/100),3857) as the_geom_webmercator,
  'DP0080005' as ethnicity,
  1 as t from tract_2010census_dp1_simple


texts: Can I save intermediary results in QGIS Processing Modeler using GRASS algorithms?

texts:
I created a 'Processing Model' in QGIS 3.0 that uses two GRASS algorithms: r.watershed and r.water.outlet. I'm interested in the output from r.water.outlet but I noticed that if I try to save intermediary results from r.watershed, the model does not work.

Ex.: In the model bellow, if I try to save the flow direction output from r.watershed, things break. If I remove the intermediary file being saved, the output from r.water.outlet works.

This question is a follow-up of the following: Second grass algorithm in QGIS modeler cannot use output from first algorithm


texts:
I would like to crop my NetCDF file for a specific location using NCO. However, I am quite lost in the documentation. So far, I've tried this:
ncks -v Precipitation_Flux -d lon,-18.751,-7.684314, -d lat,-40.9972,33.696913 C3S_concat.nc C3S_concat_cropped.nc

With error:
ncks: ERROR parsing hyperslab specification for dimension lon,-18.751,-7.684314,
Stride must be specified (and be a positive integer)

Deleting float and negative values didn't help me. And I am not understanding what is the Stride.


texts: QGIS VRT values differ from the source rasters

texts:
Background:

I have six adjacent DEMs (NED 1/3 arc-second .img files from The National Map). The min-max values among the six range from:

90.4 - 2100.0

Problem:

Using the 3.6.2 Raster > Miscellaneous > Build Virtual Raster tool, the resulting VRT's values range from:

105.7 - 1886.49

I expected the VRT's min-max value range to match the min-max range among the six rasters.

I then created a VRT from a single input raster, and again, the resulting VRT values did not match the source raster.

Varying the resample method made no difference.

Since the VRT is a "snapshot" of the input rasters, I would expect its min-max values to equal the min-max values of the its source rasters.

What am I missing?


texts:
first: my screenshot

What I want: the two fields here only with "Straße:" should be empty.
The other four results are ok.
What I have done is:

CASE WHEN  "Straße/Nr."  IS NOT NULL THEN   ('Straße: '  +  "Straße/Nr.") END

Any ideas?


texts: Buffering geometry with points in WGS84 using shapely?

texts:
I have a polygon in python shapely that I would like to buffer with a 100 foot buffer.

If the polygon is stored as a series of WGS84 lon/lat points how would I go about calculating the distance of the 100ft buffer?

I can't find anything searching google and seek guidance on the algorithm.