texts: Input data format and read multiple files into Grass 6.4.1

texts:
I am trying to read in multiple (hundreds, thousand) LiDAR ascii files into Grass.
They contain elevation data.
I think I have managed one, and I can make a map of elevations, using File, Import Raster Map, Multiple formats using GDAL.
The format is as follows:

ncols         1000
nrows         1000
xllcorner     562000.000
yllcorner     7124000.000
cellsize      1
nodata_value  -9999.0
108.93 108.92 108.97 108.99 109 109.01...
...
...

1000 lines, 1000 columns, representing 1 sq km with values at every 1 m.
Files are called

572000-7144000_dem_1m.asc
572000-7145000_dem_1m.asc
572000-7146000_dem_1m.asc

etc.

My question is, how do Y read in multiple files? I have tried concatenating them, no luck.

Thank you for your help..
regards,
Szilard


texts: QGIS 3.16 unexpected behavior when moving feature

texts:
The blue municipality boundary has a different CRS than the pink ward boundaries, and it doesn't seem possible to set them both to the same projection. If the municipality boundary CRS is changed to that of the ward boundary, it disappears.
So, the brute force solution is to simply move the blue boundary slightly west, and it will then be aligned correctly. The Vector Move Feature(s) (not Copy and Move Feature) tool is supposedly selected. Sometimes a small white circle shows up inside the Move Feature(s) icon.
It is possible to select the blue boundary, a very thin red line shows up inside it, and it's possible to move it to a new location, BUT the original blue boundary remains where it was, so now there are two blue boundaries. The select process doesn't seem to "pick up" the original blue boundary, but rather copies it and this copy is deposited in the new location.
Also, I'm not able to delete either line, so clearly my approach is fundamentally missing something.


texts: Poor performance when adding layer indexed ArcGIS Service layer to AGOL

texts:
My data sits in enterprise geoddatabases on Windows Server 2012, SQL Server manager.

I often add ArcGIS Server Web services to my AGOL maps and they perform really well. However, when I add a specified indexed layer ID when adding the service, it drastically affects the performance of the layer. e.g:

https://xxx.xxx.com/xxxx/rest/services/xxx/xxx/MapServer/0

instead of:

https://xxx.xxx.com/xxxx/rest/services/xxx/xxx/MapServer

It's confusing me because the data is coming from the same source (same enterprise geodatabase, same server, same service), but performs a lot worse when added individually, rather than as a whole service.

The reason I want to add the layer individually, is because I use it in Operations Dashboard. However, when adding the whole service to a web map, I cannot get the map within Operations Dashboard to link/filter/query with the other chart elements via their actions tab "Element does not support any events". This functionality only becomes available to me if I reference the layer ID when adding the service to a web map, or upload and use AGOL hosted feature layers instead of service layers from my enterprise geodatabase.

I created spatial indexes which did not seem to help, but i doubt indexing should make a difference, when the same source works fine when adding the service as a whole.

Could anyone tell me:
a) What could be causing this big performance drop?
b) Any recommendations to rectifying/improving this performance?
c) An alternate way to get the Operations Dashboard functionality I am after in AGOL, using enterprise-sourced ArcGIS server layers?


texts: Calculating values of field based on two others using Field Calculator in ArcGIS Desktop?

texts:
I need to combine two fields in an attribute table.  I want to do this using the field calculator.  What is the proper script to write to combine a field called "name" and "type", so that it produces one field with "name type"?

I'm using ArcGIS Desktop 10.0.


texts: Troubleshooting why ArcGIS Model works on one PC but not another?

texts:
Do you have specific expertise you could share of examples where an ArcGIS Modelbuilder has worked on one PC but not on another?

Are there any issues between different Windows versions?

There is a related question here Model not editable on a different PC

In my example I have a model that runs to completion on a Windows 8 machine. I run the models from an external hard drive. When the same model is run from the same external drive but on a different PC - Windows 7 then the model crashes.

The model can be opened, edited and verified on both Pc's.

It just crashes on one but not the other.

I have explored possible options:

ArcGIS Version and service packs - both are 10.1 SP1.
RAM - the PC on which the models work has LESS RAM than the PC on which the models crash
OS type - both machines are 64 bit.
Model build - I have re-built the model by copying and pasting into new / clean models (I did this on the Windows 8 PC though)

I am not sure that the cause is the different Windows version, it is just the main difference I can see between the two PC's

Any expertise, answers or references to sources of information would be much appreciated.

I am using ArcGIS 10.1 SP1 for Desktop Advanced.

UPDATE:

I have followed all the advice in setting up a shared Toolkit so all folder structure to Scratch, working directory, etc is set to named folders within the Toolkit. So I do not think that was the issue.

Unfortunately at the moment I do not have copies of crash report, etc to post here.

I have reported this issue to Esri, and sent them the relevant models but they were not able to reproduce the issue.

Ultimately, we upgraded our machines to 10.2.2 and we have not had the issue, although the models have been further edited a lot since that time.

I can only assume that it was some issue with 10.1.


texts:
I'm very new to QGIS, but in the past days I've created 3 different maps. I have now downloaded DEM data in ASC format, and when I open it as a raster file in a new QGIS project there is no problem, I switch it to hillshade and it displays perfectly. The problem is when I want to have this DEM in the same file as my 3 maps (polygons). It loads, but it is not visible.

If I open my 3 maps and then load my DEM, the DEM doesn't show. And if I open the DEM in a new blank file, and I load my polygons from my 3 maps, the polygons don't show. It's as if they're not compatible with each other.

I've also opened the DEM in ASC format and exported it as GeoTIFF, and afterwards I loaded the GeoTIFF format into the file where my maps are, but the the DEM still doesn't appear.


texts: 3D ground model with open source software?

texts:
I am looking for some open source software to plot some borehole logs and create a 3D ground model.

I used to work in a mining environment and used software including Datamine and Mapinfo Discover to create 3D models. I no longer have access to those programs, and am now doing a review of some GI campaign for a client and need to produce similar 3D models. I'm hoping there is an open-source software solution for this type of work.

The dataset I have is in Excel format and it consists basically of the borehole coordinates, borehole survey, geology logs and geophysical properties.


texts: How to bulk upload raster datasets to ArcGIS Online as separate discrete layers?

texts:
I have 18 raster datasets which I need to publish through ArcGIS online to allow users to turn on/off each layer and toggle transparency. I have the data stored in a file gdb as below:

I am looking at ways of doing this more efficiently than publishing 18 separate services for each layer. The service will be shared to ArcGIS Server Manager.

Can anybody help?


texts: Coastline where raster was clipped based on mask vector layer

texts:
I have clipped a DEM raster using as mask layer a (detailed) vector. Now, if I keep don't zoom in on the coastline, the two layers work very well together.

But if I zoom on any port cities, I get the pixels from the raster. Is there any technique to clip more precisely the raster, filling what is contained by the mask dropping was is outside?


texts: How unselect a feature ithat select with OpenLayers-3.8.2's selectInteraction?

texts:
I'm using OpenLayers-3.8.2. I want to select a feature and then delete it. my code is as follow:

deleteFeature = function(event){
    var feature = event.selected[0];
    myLayer.getSource().removeFeature(feature);
}

selectInteraction = new ol.interaction.Select({
    condition: ol.events.condition.singleClick,
    toggleCondition: ol.events.condition.shiftKeyOnly,
    layers: [myLayer]
})

map.addInteraction(selectInteraction);

selectInteraction.on("select", deleteFeature, this);

This code delete the feature from source but it still showed in map.

How can I unselect a feature? Or how can I redraw a layer? Or how can I delete a featrue in openlayers-3?


texts: It's seems that output of gdaltransform in meters. I would like to enter data in geography degrees and get result in degrees. Is it's possible?


texts:
I am trying to connect to an FTP site from my states GIS data center using QGIS. I am new to GIS so I'm not sure if there is a problem with the FTP, a problem with my workflow, or if this site isn't capable of being connected to.

I simply go to Layer->Add Layer->Add WMS/WMTS Layer I then create a new server connection and paste the FTP server address (ftp://www.maris.state.ms.us/ftp4/NRCS-SE_2015_DEM_IMG/)

I then try and connect and get this message.


texts: Listing datasources in MXDs using ArcPy?

texts:
I've written a tool that prompts the user for a directory path, and optionally a single mxd. It then crawls the directory for mxds (or just the single one) and writes some of the properties of the layers in each mxd to a csv file. Most of this works fine, except none of the data sources are listed.

I've seen Crawling directory and listing all datasources using ArcPy? but can't see where I am going wrong. The other properties are written just fine.

Modified the inner for loop a bit, just to see what was being found. The line:

arcpy.AddMessage('\ndataSource is: {0}'.format(lyr.dataSource))

does display the data source of all supported layers. Ran it on a larger set of mxds, with sde files, shapefiles, file gdb feature classes, web services and connections to servers, such as USDS and our own server. The only data sources that printed to the csv were for image files on the servers. The connections to these are in ArcCatalog, which is where I'm running the tool from.

import arcpy, os, fnmatch, csv

mxddirectory = arcpy.GetParameterAsText(0)
mxd_single = arcpy.GetParameterAsText(1)
outputcsvlocation = arcpy.GetParameterAsText(2)

mxd_list = []

if len(mxd_single) > 0:
    mxd_list.append(mxd_single)
else:
    for dirpath in os.walk(mxddirectory):
        for filename in dirpath[2]:
            if fnmatch.fnmatch(filename, '*.mxd'):
                mxd_list.append(os.path.join(dirpath[0], filename))

if len(mxd_list) > 0:
    outputcsv = open(outputcsvlocation, 'wb')
    writer = csv.writer(outputcsv, dialect = 'excel')

    writer.writerow(['Mxd path', 'Layer Name', 'Layer Description', 'Layer Source'])

    for mxdpath in mxd_list:
        mxdname = os.path.split(mxdpath)[1]

        try:
            mxd = arcpy.mapping.MapDocument(mxdpath)

            dfList = arcpy.mapping.ListDataFrames(mxd)

            for df in dfList:
            for lyr in arcpy.mapping.ListLayers(mxd, '', df):
                if lyr.supports('dataSource'):
                    arcpy.AddMessage('\ndataSource is: {0}'.format(lyr.dataSource))
                    layerattributes = [mxdpath.encode('utf8'), lyr.longName.encode('utf8'), lyr.description.encode('utf8'), lyr.dataSource.encode('utf8')]
                    writer.writerow(layerattributes)
                else:
                    arcpy.AddMessage('\nLayer {0} does not support dataSource.'.format(lyr.longName))

        except Exception as e:
            arcpy.AddMessage('EXCEPTION: {0}\n{1}\n{2}\n'.format(mxdpath, lyr.longName, e))

        writer.writerow(['', '', '', ''])

        del mxd

    outputcsv.close()

else:
    arcpy.AddError('No ArcMap documents found')


texts: Updating feature in ArcGIS REST service with Ajax?

texts:
I have to modifiy an app. This app must update and add data to feature in an ArcGIS Server service. I know that with ArcGIS JavaScript API is very easy but this app is build with Leaflet and others things and I can't change.

I'm trying to do it with ajax, this is my code

$.ajax({
            url: featureServiceUrl + "/updateFeatures",
            dataType: 'json',
            type: 'POST',
            data: "features:"+poiToAdd+"&f=json",
            success: function (data, textStatus, jqXHR) {
                debugger;
                //
                console.log('success');
                console.log(data);
                console.log(textStatus);
                console.log(jqXHR)
            },
            error: function (jqXHR, textStatus, errorThrown) {
                debugger;
                //
                console.log('error'); //<-- the response lands here ...
                console.log(jqXHR); //<-- the console logs the object
                console.log(textStatus); //<-- the console logs 'parsererror'
                console.log(errorThrown); //< -- the console logs 'SyntaxError {}'
            },
            complete: function (jqXHR, textStatus) {
                debugger;
                //
                console.log('complete');
                console.log(jqXHR);
                console.log(textStatus);
            }
        });

PoiToAdd is the data to modify, this data is correct because when I introduce it in this windows

the data is updated in my server, but with ajax code returns an error, that the data isn't correct. I also try with XMLHttpRequest
I have to say that this not is my code, but this work fine in other app from my company

update: function(updateFeature,callback){
    var updateparams = "features="+updateFeature+"&f=json";
    var updateurl=this.url+"/updateFeatures";
    var http;
    debugger;
    if (window.XMLHttpRequest){
        http=new XMLHttpRequest();
    }
    else {
        http=new ActiveXObject("Msxml2.XMLHTTP");
    }
    http.open("POST", updateurl, true);
    http.setRequestHeader("Content-type", "application/x-www-form-urlencoded");
    http.onreadystatechange = function() {//Call a function when the state changes.
        if(http.readyState == 4 && http.status == 200) {

            if (callback) { callback(http.responseText); }
        }
    }
    http.send(updateparams);
}, //end UPDATE

Here Updatefeatures is == to PoiToAdd and this.url is the featureservice url, now when i execute this doesn't return any error,  but the data isn't updated


texts: Updating layer style with ArcPy for ArcMap

texts:
I have written a small script.
I want to add a layer to the existing .mxd, change the style and export as pdf
My code
import arcpy
mxd = arcpy.mapping.MapDocument(r"C:\Untitled.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
arcpy.env.workspace = r"C:\new\New File Geodatabase.gdb"
arcpy.MakeFeatureLayer_management("testing", "lyr")
layer = arcpy.mapping.Layer("lyr")
arcpy.mapping.AddLayer(df, layer, "TOP")
layer2 = arcpy.mapping.Layer(r"C:\new\pointst.lyr")
arcpy.mapping.UpdateLayer(df,layer,layer2,True)
arcpy.mapping.ExportToPDF(mxd, r"C:\example.pdf", df)

The exporting is done but the style is not updated based on "pointst.lyr"(this is another layer not currently in the mxd).
I have used arcpy.mapping.UpdateLayer and "true" since I needed to just update the symbology.


texts: How to specify default value for Rotation in SLD?

texts:
Am using geotools to render a set of point features. I have SLD for styling point features. The SLD uses an expression to calculate Rotation value based on Angle attribute value of feature. The features that have Angle value NULL are not rendered. I would like to specify default angle value of 0 when Angle attribute value is NULL.

Here is part of current SLD:

<sld:Rule>
<sld:Name>rule1</sld:Name>
<sld:PointSymbolizer>
<sld:Graphic>
<sld:ExternalGraphic>
<sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:href="mysymbol.svg" xlink:type="simple"/>
<sld:Format>image/svg</sld:Format>
</sld:ExternalGraphic>
<sld:Opacity>
<ogc:Literal>1.0</ogc:Literal>
</sld:Opacity>
<sld:Size>
<ogc:Literal>11.0</ogc:Literal>
</sld:Size>
<sld:Rotation>
<ogc:Sub>
<ogc:Literal>360.0</ogc:Literal>
<ogc:PropertyName>Angle</ogc:PropertyName>
</ogc:Sub>
</sld:Rotation>
</sld:Graphic>
</sld:PointSymbolizer>
</sld:Rule>


texts: Oracle to PostgreSQL migration with ogr2ogr only imports first geoloc column

texts:
I transfer Oracle tables to PostGIS using Oracle Wallet and

 ogr2ogr -f "PostgreSQL" PG:"host=<host> port=5432 user=<user> dbname=<dbname> password=<pwd>" OCI:<USERNAME>/<PASSWORD>@<SID>:<table>

It transfers all regular columns and the first geoloc column. Tables with multiple geoloc columns are not transferred correctly. I can create pg tables with multiple geoloc columns by joining several tables within Pg client.

Is there a way to copy multiple geoloc columns with ogr2ogr?


texts:
Background, in case it helps: My objective is to (eventually) fully automate the mass digitization of grain silos. To this end, I have set out to first find a way to automatically identify silos in a given raster image.

Question: Is there any way to automate the generation of training samples for the purpose of supervised classification in Arc, QGIS, or any other program?

I'm pawing through the literature and mainly finding things about convolutional neural networks and machine learning which, frankly, is above my pay grade. It's sadly looking like I could potentially be better off just with unsupervised classification, at the cost of accuracy and specificity.

As always, happy for any input or literature review recommendations!


texts: Removing label arrows using Easy Custom Labeling tool?

texts:
I am labeling a lot of features and the Easy Custom Labeling tool is helping me to move every label to the desired position. I also like that the tool adds arrows for labels that I moved but sometimes I am just moving them some few pixels and do not want an arrow to be created.

I am now wondering how I can get rid of a selection of arrows.
For now, I can only think of creating two layers: one for labels without arrows and one for the ones with arrows.


texts: Change of Internal polygon area after intersecting polygon


texts: How do I find/define a projection that gives me exact relative positions inside an area?

texts:
Assuming that a set of positions are all inside a 5km square area:
How do I find/define a projection to reproject the positions, so that the relative positions between the now reprojected positions are exact?

This would mean that for any pair of objects inside that 5km square area; If I do distance = obj1.EastingNorthing - obj2.EastingNorthing I will get a distance vector whose value has less than 1cm of error, in relation to the true distance.

Is there a single projection that I can use globally?
Does it have to be a different projection for each square? How do I build it?


texts: Exporting entire ImageCollection in Google Earth Engine using geetools for JavaScript or Python for Windows?

texts:
I am very new to ee and am trying to export a collection to my drive, after masking for clouds and clipping to myregion... tried to follow the advice in this answer https://stackoverflow.com/questions/35991403/pip-install-unroll-python-setup-py-egg-info-failed-with-error-code-1 where the col2drive function is imported from geetools. But I seem not to have access to the repo, it doesn't show under my read module. Took me a while to find out because the error said is not a function, and thought I was not creating the object correctly.
So I copied the hardcoded version posted here Export each image from a collection in Google Earth Engine to my own repo and tried to load it from there, did not work. I get again and error saying col2drive is not a function... when I do  var exportfunction = tools.col2drive; and  print(exportfuncion)... it prints "undefined". So is not loading it, I think.

Then I found this other option Google Earth Engine: Export an entire collection, running it directly from python API. Installed python in my computer (Windows10, 64 bits), and when trying to install geetools it needs a module  'fcntl', which apparently does not exist for Windows... there is a suggestion of using  win32api, but then again is for windows32 and the advice says that it is important to understand what  fcntl is doing in each case before substituing... In any case I have no idea where can I substitue it.
Could someone help, I did try.

Here is my (stitched) code in case is useful

var geometry = /* color: #98ff00 */geometry;

//Add tools to export collection
// load tools from Fito Principe
//var tools = require('users/fitoprincipe/geetools:tools');
var tools = require('users/sulazul/Rep:col2Drive')
// add col2drive function
var exportfunction = tools.col2Drive;
print(exportfunction)

// Load Landsat 8 imagery and filter it
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                   .filter(ee.Filter.eq('WRS_PATH', 33))
                   .filter(ee.Filter.eq('WRS_ROW', 41))
                   .set('SENSOR_ID', 'OLI_TIRS')
                   .sort("CLOUD_COVER", false);
var count = collection.size()
print("Coleccion", count)

//some more code here to maskclouds and clip to myregion

// Map the function .
var filtered_1 = filtered.map(maskL8sr);
// export collection to google drive
exportfunction(filtered_1,"Landsat8", {
  scale: 30,
  region: geometry.getInfo()["coordinates"],
  type: 'float'
});


texts: Clipping line layer to every feature of polygon layer in QGIS?

texts: I have a line layer (roads) and a polygon layer which represents districts/counties. I want to clip the roads according to every one different county. That means that I must repeat the procedure sixty times as the total count of districts is 60. (I must select one different district (feature of polygon) every time and then Vector -> Geoprocessing -> Clip). Is there any possible way this procedure to become automatic? I work in QGIS (and ArcGIS).


texts: Creating second map without labels using QGIS?

texts:
I want to remove the labels in the second, smaller map. So far when I switch off all the labels and refresh only the smaller map, of course, the labels are not shown. But it seems like QGIS updates also the other map when I save the map as pdf.

Is there a way to show in this map only the labels from the Rettungsleitpunkte and exclude all the others?


texts: Where Can I find US Roads? “OverWeight Trucking Business”

texts:
I tried OSM, Google Maps, TIGER Database... But there's nothing like how the permits of the OverWeight Trucks roads are issued.

Any idea about this?


texts: I have a list of .tif files which have the wrong projection system (NAD_1983_UTM_Zone_13N). However, I want to have those files in (WGS_1984). I want to do this by a script that can transform tif files from one projection system to another. I have found many scripts online, but none seem to work for me. anyone a script or idea?


texts: Openlayers: with a list of features, how can I organize them by distance from a coordinate

texts:
Let's say I'm using openlayers (5) and I have an array of features, they can all be from different layer sources, and I have a set of coordinates.  What is the best way to sort that array so that the closest feature to that coordinate is at the top, furthest is at the bottom?

I was looking at https://openlayers.org/en/latest/apidoc/module-ol_source_Cluster-Cluster.html#getClosestFeatureToCoordinate originally, but that only gives me 1 feature and not really what I'm looking for.


texts: How to compute population density of certain type with Arcmap table

texts:
I have to find areas on a map where the yearly income is less than a certain threshold, I have all of the data laid out in a table but I am not sure how to process it.

The density for a statistical area can be defined as the ratio of:
 - (the number people who earn less than x in area1) / (the area of the
   area1 (in hectare))

Statistical areas are classified as ‘high density’ if they have a density greater than 3 people who earn x or less per hectare.

table is laid out like this: $x to $y amount = 20 people, for example. (so they have all been categorised already)

I have worked out how to create a list of the density by simply doing the division above, but how would I classify these numbers into < 3?

Also, how would I go about ignoring, or not displaying the areas that do not suit my needs (eg. > 3)


texts: Using ExtractData in ArcGIS API for Javascript?

texts:
I cannot find any sample using ExtractData API.

I am trying to extract data from hosted feature layers using ExtractData JS API (not widget). The job that I submit fail with  the following message.

  Error: [{"type":"esriJobMessageTypeError","description":"ERROR 000735:
  inputLayers: Value is
  required"},{"type":"esriJobMessageTypeError","description":"Failed."}]

Following is part of the request.  The 'instance' object has root url and credentials. All layers are accessible and also I had made them public for testing purpose. I have spend lot of time trying to resolve this issue. There is no other meaningful information coming out of GP Service.

const job = await instance.get('/submitJob', {
        params: {
          inputLayers: [
             { url: soilURL, serviceToken: user.creds.token },
             { url: forestURL, serviceToken: user.creds.token },
             { url: homesiteURL, serviceToken: user.creds.token},
          ],

          extent: {
            url: parcelURL,
            serviceToken: user.creds.token,
            filter: parcelFilter,
          },
          clip: true,
          dataFormat: 'FILEGEODATABASE ',
          outSR: 102100,
          outputName: {
            title: 'ACT 319 clip data',
            tag: 'clip data',
            snippet: 'clip data',
            description:
              'Clip and download soil type, home site and forest layers.',
          },
        },
      });


texts: Losing data when updating ArcGIS for Desktop?

texts:
I would like to update my ArcGIS from version 10.3 to 10.4 (in particular, I am interested in ArcMap and arcpy).

I read that on the ESRI webpage referring to the update from 10.1 to 10.2 that any customizations, add-ons, or third-party applications that were installed for ArcGIS 10.1 must be reapplied for ArcGIS 10.2.

If I update, my ArcGIS folder (the folder created during ArcGIS installation, C:\Users\Myusername\Documents\ArcGIS, where I stored all the shapefiles, layers, bookmarks etc.) will be deleted? Or is there any chance I will lose my data? If yes, how can I prevent this?


texts: Leaflet.VectorGrid throws error Input data is not a valid GeoJSON object

texts:
This is a cross-post from SO. I am trying to use Leaflet.VectorGrid.Slicer, but I am getting error:
Input data is not a valid GeoJSON object

codesandbox
VectorGrid.js
import { useEffect } from "react";
import { useLeafletContext } from "@react-leaflet/core";
import L from "leaflet";
import geojsonvt from "geojson-vt";
import "leaflet.vectorgrid";
import { featureCollection, feature } from "@turf/helpers";

export default function VectorGrid({ assets }) {
  const { layerContainer, map } = useLeafletContext();

  const geoJsons = assets
    .filter((asset) => asset.mapPolygon)
    .map((asset) => feature(asset.mapPolygon.geometry));

  const slicerInput = featureCollection(geoJsons);

  console.log(slicerInput);

  const data = geojsonvt(slicerInput);

  const vectorGrid = L.vectorGrid.slicer(data, {
    rendererFactory: L.canvas.tile,
    vectorTileLayerStyles: {
      sliced: {
        fillColor: "transparent",
        color: "blue",
        weight: 0.5,
      },
    },
    maxZoom: 22,
    indexMaxZoom: 5, // max zoom in the initial tile index
    interactive: true,
  });

  const container = layerContainer || map;

  useEffect(() => {
    container.addLayer(vectorGrid);
    return () => {
      container.removeLayer(vectorGrid);
    };
  }, [container, vectorGrid]);

  return null;
}

temp_data.js
export const dat = [
  {
    assetid: 1,
    id: 1,
    mapPolygon: {
      type: "Feature",
      geometry: {
        type: "Polygon",
        coordinates: [
          [
            [-104.98948, 39.64244],
            [-104.98953, 39.64244],
            [-104.98953, 39.64236],
            [-104.98918, 39.64236],
            [-104.98918, 39.64247],
            [-104.98948, 39.64247],
            [-104.98948, 39.64244],
          ],
        ],
      },
    },
    latitude: null,
    longitude: null,
  }, {
    assetid: 3,
    id: 3,
    mapPolygon: {
      type: "Feature",
      geometry: {
        type: "Polygon",
        coordinates: [
          [
            [-104.98793, 39.6444],
            [-104.98816, 39.64439],
            [-104.98816, 39.64438],
            [-104.98825, 39.64437],
            [-104.98825, 39.64432],
            [-104.98826, 39.64431],
            [-104.98826, 39.64429],
            [-104.98828, 39.64429],
            [-104.98827, 39.64422],
            [-104.988, 39.64422],
            [-104.98799, 39.64429],
            [-104.98793, 39.64429],
            [-104.98793, 39.6444],
          ],
        ],
      },
    },
    latitude: null,
    longitude: null,
  }]


texts: Exported PDF resulting in misaligned Bing Aerial and vector layer using QGIS?

texts:
I'm working with QGIS 2.8.1. I have a project with CRS: 3947 (Central-Western France) and several layers in the same CRS. I have not activated the 'On the fly' CRS transformation.

I add the Bing Aerial layer through OpenLayers plugin. My project then automatically goes in CRS 3857 (OTF), despite the settings explained before... This is not a big issue because my vector layers display automatically follows the Bing Aerial and align perfectly.

I then create a new composer with a map window. The vector layers are perfectly aligned on the Bing Aerial background. The issue comes when I export the composer as a PDF or PNG, resulting in a big misalignment between vector layers and Bing Aerial... I don't understand why the printing process modifies what is displayed in the composer.

Any idea about this issue and how I could solve it?


texts:
I am using leaflet to build a layer control with three layers all of them shouldn't be on the map simultaneously. A layer should only appear on the map with a click. the code below doesn't work

   var attractionnyc = L.layerGroup();
   var museumsnyc = L.layerGroup();
   var hotelsnyc = L.layerGroup();

    L.control.layers ({}, {
         'Attraction': attractionnyc,
    },  {
        collapsed: False,
    }).addTo(map);

      L.control.layers ({}, {
         'Museum': museumsnyc,
    },  {
        collapsed: False,
    }).addTo(map);

        L.control.layers ({}, {
         'Hotels': hotelsnyc,
    },  {
        collapsed: False,
    }).addTo(map);


texts: How to store attributes of a large shapefile in a list?

texts:
I am trying to store the attributes of a shapefile into a list. With following code it works for smaller shapes:

addA = []
While provider.nextFeature(feat):
    attrs = feat.attributeMap()
    for (k,attr) in attrs.iteritems():
        addA.append(attr.toString())

If I try to store a large shape my pc runs out of memory or freezes. I tried to put the attributes in a blist and in a numpy array, but nothing really worked (or probably I did not understand how to set up a proper numpy array...). Is there a solution for this problem?


texts: Counting number of pixels in each tile in Earth Engine

texts:
I have 0.02° tiles across some regions of interest in Africa. I am trying to find the percent of 300m ESA land cover pixels that are classified as cropland (classification 10-30) for each tile, as well as the total number of ESA non-NA pixels. I thought one way to do this would be to count the within-tile pixels classified as cropland and then count the total number of ESA pixels within each tile. I did this as follows:
//Select only the 2001 ESA image
var lulc_01 = lulcYrSub.select(["lc_2001"]);

//Reclassify the ESA image such that all pixels with classification 10-30 are 1s and the rest zeroes
var ag_2001 = lulc_01.gte(10).and(lulc_01.lte(30));

//Reclassify ESA image so that all non-NA pixels are 1's
var tot_2001 = lulc_01.gte(10).and(lulc_01.lte(220));

//Function adds up total number of cropland pixels within each tile
var nAg_2001 = tiles.map(function(feature) {
  return feature.set(ag_2001.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: feature.geometry(),
    scale: 30
  }));
});

//Function adds up total number of pixels within each tile
var nPix_2001 = nAg_2001.map(function(feature) {
  return feature.set(tot_2001.reduceRegion({
    reducer: ee.Reducer.sum(),
    geometry: feature.geometry(),
    scale: 30
  }));
});

But the resulting pixel counts seem impossible given the tile size. My tiles are about 2.22 km by 2.22 km. There should be approximately 50-60 300-meter pixels within each tile. Yet one of my tiles allegedly has 1,900 pixels in it! What is going wrong? Can you recommend an alternative way to get the percent crop cover as well as the total number of pixels within each tile?  Full code here.


texts: Preserving links when using ArcMap and going between Different Computers?

texts:
I am a coop student who is close to finishing my second GIS coop term, and once again coming across the issue the fact that when my computer is wiped all the links to my files in ArcMap will be broken. For example, I created a few different ModelBuilder processes on my computer, if I put it on the shared local drive I will have to go through all the input and output files in the model and rename and remake folders.

Is there a way to make the files pathways dynamic, or some way of saving all the files to a process so that I will not have to change all the pathways on every different drive?


texts: Converting ee.Date object to client string in Google Earth Engine?

texts:
I'm trying to export all the images in an ImageCollection to Google Drive and by default name the output image with the date using format() function. However, when I try to name a export task by a date, the name in tasks is something like:

ee.String({ "type": "Invocation", "arguments": { "date": { "type": "Invocation", "arguments": { "value": { "type": "Invocation", "arguments": { "object": { "type": "ArgumentRef", "value": "_MAPPING_VAR_0_0" }, "property": "system:time_start" }, "functionName": "Element.get" } }, "functionName": "Date" }, "format": "YYYY-MM-DD" }, "functionName": "Date.format" })

Note the date values are gone and only some of the images are exported. Here's the link to the code

All I need is only a formatted string from a date object.

Is there a simpler way to export an image collection?


texts: Setting GeoServer data directory location?

texts:
I want to run GeoServer 2.14 with Jetty on a Ubuntu 18.04.
I copied the war-file in the webapps folder and GeoServer is responding.

Unfortunately I can not set the data directory location as described in the GeoServer docs.

I appended the following line to /etc/environment:

GEOSERVER_DATA_DIR="/opt/geoserver_data"

After logging in to the system the variable is set as expected:

#echo $GEOSERVER_DATA_DIR
/opt/geoserver_data

I also created the folder /opt/geoserver_data:

#ls -la /opt/geoserver_data
total 8
drwxr-xr-x 2 root root 4096 Nov 12 15:20 .
drwxr-xr-x 5 root root 4096 Nov 12 15:20 ..

When I restart Jetty, the first two GeoServer log lines are:

INFO [geoserver.platform] - Falling back to embedded data directory: /tmp/jetty-0.0.0.0-8080-geoserver.war-_geoserver-any-6077051313967886521.dir/webapp/data
INFO [geoserver.platform] - Falling back to embedded data directory: /tmp/jetty-0.0.0.0-8080-geoserver.war-_geoserver-any-6077051313967886521.dir/webapp/data

This is strange because according to the GeoServer source code, there should be a warning that the given folder does not exist or is not writable.

Obviously I am doing something wrong, but I have no idea what?

Based on Ian Turton's answer, I ensured that jetty is not run as root.
Furthermore I start jetty now with -DGEOSERVER_DATA_DIR=/opt/geoserver_data, which is easier to handle than a system variable.

At the moment I still get the following warning messages:

WARN [geoserver.platform] - Found Java environment variable GEOSERVER_DATA_DIR set to /opt/geoserver_data , which is not writeable
INFO [geoserver.platform] - Falling back to embedded data directory: /tmp/jetty-0.0.0.0-8080-geoserver.war-_geoserver-any-1540535100624462590.dir/webapp/data

But at least the expected warning is shown. Although the fallback data directory still looks odd.


texts: How to add layer/raster with specific position in the QGIS layer list?

texts:
I use addRasterLayer(rast,'Name') but this adds the raster to the top.
It is possible to set position when I add raster at once? Is it possible??


texts: Add custom data into Nominatim geocoding search?

texts:
I have installed Nominatim on my local server to receive addresses from long/lat. I would like to add custom data to the search itself. I would like to store locations of the rooms of our buildings.  Is this possible?

We are storing polygon data for different types (layers) of objects. How can those be returned by Nominatim?
E.x.:
 - return data from layer: water_protection_areas
 - return room polygons for specific address for 4th floor
 - return room id(s) for specific geocode

(by Mihail-Klr, we are researching on this topic, if we should find a solution, I'll post it)


texts: PostGIS Shapefile DBF Loader 2.2 in pgAdmin III greyed out - why?

texts:
I have two issues am I am wondering if they are related.

1: I can't get plugins to enable in pgAdmin III
I have two PostgreSQL 9.5 instances. In one, the pgAdmin III plugin "PostGIS Shapefile DBF Loader 2.2" works when I highlight my database. In the other instance that plugin menu choice is greyed out. I've been comparing the two instances and I can't figure out what is different. In both cases I am connected to Postgres as a user who is not the owner of the database. PostGIS is installed in both. For the plugin, what are the requirements for it to be active (not greyed out)?

Additional info 3/28/2017:
Note: my PostgreSQL is on Linux and my pgAdmin III is on Windows 7.
I can successfully run the shapefile loader shp2pgsql-gui.exe found in Program Files> PostgreSQL > 9.5 > bin > postgisgui. I can verify that the shapefile has loaded well by viewing it in QGIS.  The PostGIS Shapefile DBF Loader 2.2 and the PSQL Console are both greyed out when I select the database. I am connecting as afreeman which is a superuser.  I also tried connecting as enterprisedb which is a superuser and owner of the public schema.
My "PG bin path" as seen in Options/Browser/Binary paths is C:/Program Files/PostgreSQL/9.5/bin

2: I can't see a listing of geospatial layers in QGIS
My QGIS will connect to one PostGIS instance but not list layers in any schemas. I was able to load shapefiles using the QGIS DBManager into schemas afreeman or public and I can drag those layers into QGIS map and they draw OK. However, choose Layer menu, Add Layer, Add PostGIS layer and connecting to the PostGIS database doesn't list anything. I can do this successfully with a different PostGIS instance. Same version (9.5). One difference though as that in the successful db the public schema is owned by postgres. In the unsuccessful one the public schema is owned by enterprisedb.

Does it matter?

Why would a test by successful but no layers appear?

Additional info 3/28/2017:
show search_path; gives me: "$user", public, topology

afreeman is a superuser. I'm running PostGIS 2.1 However, I've also tried connecting as enterprisedb which is the owner of the public schema, same issue with not seeing layers. When I connect as enterprisedb or afreeman and test the connection that is successful. Just no layers in the list. Ran the SQL offered below but no effect.  PostGIS is installed in the public schema.


texts: Iterating a model to create 3 buffers around three points and find their intersection (trilateration) in arcgis 10

texts:
This past summer I recorded the location of a large number of plants I was working with by calculating the distance from each plant to each of 3 tripods for which I have GPS coordinates. I'm trying to use arcgis 10 to reconstruct the location of each of these plants. I don't have the bearing from the focal plant to the reference points, only the distance. So as I understand it, this is a trilateration problem.

I've created the following model:

So I have the three GPS reference points as separate features, and each time I enter the distance from the focal plant to each of these reference points. It creates the three buffers, finds their intersection, and plots a point at the center of the intersection polygon.

Although this works, it is incredibly time consuming (~5 hours to get the first of four populations). I have around 1,300 more samples in three separate populations that I need to map in this way.

I cannot figure out how to make an iterative version of this model, however, in large part due to my inexperience with gis. I have a hard time searching the help documents because the jargon is unfamiliar to me, so I think I just lack the correct keywords.

I have seen this post, but as I searched to see how to convert between UTM and the coordinate system used there, it seemed that would be potentially problematic, and I feel like I'm 95% of the way to the solution with the model I already have.

The most obvious way this should work to me is to have an input feature that contains all of the reference and buffer distance data, and iterating through sets of three, or through a set of files each of which contains three points. However, when I've experimented with entering a table like this:

the intersection tool no longer outputs a single polygon where all three buffers overlap, but instead multiple polygons where there is an overlap between at least 2 of the buffers, like this:

I assume that this is because the three buffers are in a single feature instead of three different features, but why this would create a difference in the behavior of the intersection tool passes my understanding.

Furthermore, I don't seem to understand how to make iterators operate properly--my attempts at using the row or feature iterators always result in only one output representing the last row, even when I use dynamic naming of the output.

Is there some way to create an iterative version of this model, or should I just give up on it and go back to square one? I've struck out asking all of the GIS experts on our campus, none of whom seem to have much experience building models.

Thank you! In reading back through this I appreciate that it is very convoluted. I apologize for the length and for my ignorance--this is what happens when a geneticist pretends to be an ecologist.


texts:
I use pyqgis to create standalone GIS apps.

In my system I have only pythonpath from QGIS C:\OSGeo4W64\bin\python.exe I don't have any other python version.
in my installation folder I have a C:\OSGeo4W64\apps\Python27 and C:\OSGeo4W64\apps\Python27\Lib\site-packages folders.

can I install packages in this folders using pip ?

for example If I have install standard python and from cmd use pip then I can easy packages on this folder C:\Python27.


texts: Using feature class name as name in output (buffer_analysis)

texts:
I'm pretty new in python, and I'm trying to iterate through a list of feature classes, and perform a buffer on each feature class. So I've done this:

import arcpy
from arcpy import env
env.workspace = "C:/EsriPress/Python/Data/Midterm/Oil_project.gdb"
env.overwriteOutput = True
fclist = arcpy.ListFeatureClasses("", "Polyline", "Well_Data")

for fc in fclist
buffer = arcpy.Buffer_analysis(fc, fcname_buffer, "75 FEET")

However, I want each output file name to be the original feature class name with _buffer at the end. For example, the name of the first feature class in the list is Well_A_BB05.shp. I want the output to be "Well_A_BB05_buffer". I can't get this to work.


texts: What exactly does "enable geodatabase" perform on a SQL database

texts: This is for those extra technical folks out there. What does "enable geodatabase" do to the structure and capabilities of a Microsoft 2012 SQL database. Since I'm a GIS guy and not a SQL guy, I need the technical capabilities and enhanced functions enumerated for my SQL DBAs. Could someone provide some technical input?


texts: How can I clip the NDVI TIFF file using rasterio mask and fillnodata. I tried using the following approach but edges are not getting smooth

texts:
While clipping TIFF file edges are not getting clipped and smooth and cropped.
with fiona.open(TiffFilePath, "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]

 with rasterio.open(TiffFilePath) as src:
                out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
                out_image = rasterio.fill.fillnodata(image=out_image, mask=out_image, max_search_distance=300.0, smoothing_iterations=0)
                out_meta = src.meta.copy()
                out_meta.update({"driver": "GTiff",
                             "height": out_image.shape[1],
                             "width": out_image.shape[2],
                             "transform": out_transform})
     with rasterio.open(ClippedTiffFilePath, "w", **out_meta) as dest:
                dest.write(out_image)

**Converting tiff from epsg:4326 with sentinalsat2 satellite**

try:
            ### Change EPSG values to 4326
            ###File generated after running calculateNDVI.py
            ### NDVIData.tiff is file generated after running calculateNDVI.py
            with rasterio.open(TiffGenerationFile) as src:
                transform, width, height = calculate_default_transform(
                    src.crs, dst_crs, src.width, src.height, *src.bounds)
                kwargs = src.meta.copy()
                kwargs.update({
                    'crs': dst_crs,
                    'transform': transform,
                    'width': width,
                    'height': height
                })

                ###Resulting file after changing EPSG values
                with rasterio.open(TiffGenerationFile, 'w', **kwargs) as dst:
                    for i in range(1, src.count + 1):
                        reproject(
                            source=rasterio.band(src, i),
                            destination=rasterio.band(dst, i),
                            src_transform=src.transform,
                            src_crs=src.crs,
                            dst_transform=transform,
                            dst_crs=dst_crs,
                            resampling=Resampling.q3)
            return True


texts: How to select a part of a land cover type a certain amount of feet away from another land cover type

texts:
I have raster Landcover type Urban combined with Forest all within the county raster perimeter. I want to make a layer of only the forest cells that are more than 500 feet from an urban cell, so that all forest cells that share space with urban less than 500 ft are removed.I suppose I need to create a grid of distance from urban land, but I don't know how to input "from" into a raster calculator operation.

I want the forest/shrub that is more than 500 feet from urban land


texts: Export all QGIS-project-LayerDefinitions (qlr) using python

texts:
I try to export all QGIS-LayerTreeGroups and LayerTreeLayer-definitions (as *.qlr) using following python code (recursive search to search the whole layerTree; the path to save the *.qlr is build from group- and layer-names).

I do not know how to use the "exportLayerDefinition(..,..)".
Python says "[...] AttributeError: 'QgsLayerTreeLayer' object has no attribute 'exportLayerDefinition'" (line 8 +13)

(Using QGIS 2.14 / 2.16)
My Code Sample is

def deeperExport(laTrGr, treeName):
    pos = 0
    for child in laTrGr.children():
        if isinstance(child, QgsLayerTreeGroup):
            treeName += "/" + child.name()
            deeperExport(child,treeName)
            path= treeName + ".qlr"
            child.exportLayerDefinition(path, child.children())
            pos = treeName.rfind("/")
            treeName = treeName[0:pos]
        elif isinstance(child, QgsLayerTreeLayer):
            path = treeName+"/" + child.layerName()+".qlr"
            child.exportLayerDefinition(path, child.children())

#main
root = QgsProject.instance().layerTreeRoot()
print root
treeName = "c:/..."
deeperExport(root, treeName)


texts: How to change the edit widget for fields in a layer with PyQGIS?

texts:
Is possible change the edit widget with PyQGIS?

I have made an extension that loads a PostGIS layer and I want to edit the date fields with the date widget by default without any change in the settings of layer ( "layer properties dialog").


texts: Earth Engine: Map over list versus ee.list?

texts:
I want to use a list (of image names) as input to map, and obtain a list of images. Curiously, this works when the input is a (JavaScript) list, but not when it is an ee.List()!? In other words:

ee.List(inputs).map(function(x) {return ee.Image(x)}): does not work (returns list of strings)
inputs.map(function(x) {return ee.Image(x)}): works (returns list of images)!

Why? What is wrong in doing ee.List() first? What is ee.List().map() doing compared to list.map()? The problem is that I want to use the ee.List() approach as it seems Python has no native list.map() function.

Pseudo-example (my case involves two distinct images not available elsewhere in an ImageCollection):

var list_years = ['USDA/NASS/CDL/2016', 'USDA/NASS/CDL/2017']

var CDLs = ee.List(list_years).map(function(x) { return ee.Image(x)})
var CDLs2 = list_years.map(function(x) { return ee.Image(x)})

print(CDLs)
print(CDLs2, "CDLs2")

Link to code: https://code.earthengine.google.com/7425530dc7eadd49ab87c2dfff416d0a


texts: Why ArcGIS application does not have label manager?

texts: Why is there no label manager in the pictured application of ArcGIS for Desktop?


texts: Difference between Spatial Join Tool and Spatial Join Function in ArcGIS Desktop?

texts:
I am curious to know if anyone can shed some light on the difference between ArcGIS' Spatial Join geoprocessing tool (found in Analysis -> Overlay toolbox) and the Spatial Join 'wizard' (found by right-clicking a layer -> Join -> then selecting "Join data from another layer based on spatial location" from the first drop-down).

There have been multiple times where the tool has worked and the other has failed and vice versa, especially when using ModelBuilder and Python.


texts: Can I make QGIS export one thing at the time?

texts:
While I am exporting a large amounts of TIFF files for my work, I am having a small bottleneck in my workflow. Up until now I have devided the TIFF images in smaller pieces so that my computer can process them (a project can contain 16 GeoTIFF images or more). The exports go well one by one, however to keep productive on a larger scale during the day. I like to export all the pieces at once, and put them in a the waiting list to export. This means I can focus on other work for a few hours and come back when the pc finishes running. However at the moment there are multiple files being exported at the time. Example 5 files being exported and the rest of the 16 are in the queue.
However I notice that my PC runs hot, and finds it hard to export more then one file at the time.
Is there a way to change the waitnglist settings in QGIS to that only one file is processed at the time, and the rest is in the waiting list?
Technical details:
Using version 3.4 of QGIS, on a Windows pc.
Also I am not a coder.


texts: How to add multiple markers on map and popup different information using OpenLayers?

texts: Can anybody give me some examples of adding multiple markers on the map?


texts: Get feature value to use in print composer label

texts:
I'm trying to get the value of a feature (the area) to use as label in the print composer, but I'm not getting much success.
This is the expression I'm using in the text label in composer:
attribute(get_feature_by_id ('POLIGONAL_FNHIS2008A_8babdf78_50dc_480f_b811_0da927b7c9a5',1),'area')

The program shows me this message: No feature was found on this layer to evaluate the expression.
But the image bellow contains the information about my layer I'm referring in the expression, showing the fields and the values. What am I doing wrong?


texts: Converting CSV file with local coordinates for project area to UTM zone 29N using QGIS

texts:
I have a csv file with local coordinates for a project area and want to convert to UTM zone 29N.
How do I do it with QGIS?
The files plot outside the project area, when imported into QGIS, (X-5022.867, Y-7259.176, Z-242.878) CRS is EPSG
:32629.


texts:
We're currently looking for a Land Use dataset for part of England (UK) and were wondering what datasets are out there, either free or commercial.

One's we've found so far:

http://www.infoterra.co.uk/landbase
http://www.geoinformationgroup.co.uk/products/national-lu
http://data.gov.uk/dataset/land_use_statistics_generalised_land_use_database (more generic statistics)

Important: We're not looking to generate this data ourselves, We're only interested in nicely pre-packaged datasets.


texts: Grid registration issue converting from GeoTIFF to Esri ASCII Raster format?

texts:
I've got a 3 arc second digital terrain model represented as a 32-bit GeoTIFF that looks like this:

$ gdal_info gom03_v1_0.tif

Driver: GTiff/GeoTIFF
Files: gom03_v1_0.tif
Size is 10200, 7800
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-71.499580383300781,45.999583333330001)
Pixel Size = (0.000833333333333,-0.000833333333332)
Metadata:
  TIFFTAG_SOFTWARE=nc2geotiff
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -71.4995804,  45.9995833) ( 71d29'58.49"W, 45d59'58.50"N)
Lower Left  ( -71.4995804,  39.4995833) ( 71d29'58.49"W, 39d29'58.50"N)
Upper Right ( -62.9995804,  45.9995833) ( 62d59'58.49"W, 45d59'58.50"N)
Lower Right ( -62.9995804,  39.4995833) ( 62d59'58.49"W, 39d29'58.50"N)
Center      ( -67.2495804,  42.7495833) ( 67d14'58.49"W, 42d44'58.50"N)
Band 1 Block=10200x1 Type=Float32, ColorInterp=Gray

Note that the AREA_OR_POINT attribute is set to AREA.

When I convert this to ARC ASCII Grid format using GDAL:

gdal_translate gom03_v1_0.tif -of AAIGrid test.asc

and look at the top of the resulting ASCII file, I see:

$ more test.asc

ncols        10200
nrows        7800
xllcorner    -71.499580383301
yllcorner    39.499583333335
cellsize     0.000833333333
 315.814208984375 314.712371826171875 ...

These xllcorner,yllcorner values are the same as the GeoTIFF Lower Left values above.  From the official description of the ESRI ASCII Raster Format, however, either xllcorner,yllcorner or xllcenter,yllcenter can be specified.  Since the GeoTIFF values represented AREA (not POINT), isn't the specification of xllcorner and yllcorner incorrect?

Shouldn't the ASCII file in this case really read as below?

$ more test.asc

ncols        10200
nrows        7800
xllcenter    -71.499580383301
yllcenter    39.499583333335
cellsize     0.000833333333
 315.814208984375 314.712371826171875 ...


texts:
How do you decide datum size (axies a & b) and datum position relative to WGS84 (dx,dy,dz,rx,ry,rz) for a country if you were to create a new datum best fit for the country?
Example, is datum surface needed to be below the lowest elevation always or can be placed at average elevation of the country topography..etc?


texts: How to import Excel data and display it on a state map in QGIS?

texts:
I want to import summarized data from Excel into QGIS 2.4. The date consist of US states and Canadian provinces paired up with a percentage. The data looks like this:

State     %
CO       50.4
IN       34.1
KY       10.1
OH       5.2
TN       24.8

How do I import this data into QGIS to display it on a map? I would like the state to be a solid color that varies with the percentage.
I would like to either have % transparency or fade from one color to the other. Also I am starting with the Natural Earth Quickstarter Kit, so I have state borders on my map.

If there are tutorials on this, please share.


texts: How to configure Saga Tools on Qgis script

texts:
I'm trying to use the process saga:paramemeltonruggednessnumber as following :
melton_saga = processing.run("saga:parememeltonruggednessnumber",
         {'DEM' : path_filled,
          'AREA' : path_area,
          'ZMAX' : path_zmax,
          'MRN' : melton_layer})

It makes me the following error :
call saga_cmd ta_hydrology "ParameMelton Ruggedness Number" -DEM "C:\Users\UTILIS~1\AppData\Local\Temp\tmpff7k_q4j\filled.sgrd" -AREA "C:\Users\UTILIS~1\AppData\Local\Temp\tmpff7k_q4j\area.tif" -ZMAX "C:\Users\UTILIS~1\AppData\Local\Temp\tmpff7k_q4j\zmax.tif" -MRN "C:/Users/Utilisateur/AppData/Local/Temp/processing_HIbIsN/94cd52427f7048f5b877831b3172d028/melton.tif"
             Error: select a tool

The Saga Gis command line must be like :
Usage: saga_cmd ta_hydrology 23 [-DEM <str>] [-AREA <str>] [-ZMAX <str>] [-MRN <str>]
  -DEM:<str>    Elevation
    Grid (input)
  -AREA:<str>   Catchment Area
    Grid (output)
  -ZMAX:<str>   Maximum Height
    Grid (output)
  -MRN:<str>    Melton Ruggedness Number
    Grid (output)

The problem here is that i'm missing the number 23, how can I implant it ?
Rgds,


texts: Invalid references created by pgr_createTopology

texts:
I am trying pgr_createTopology on a simple network and it is leaving invalid node references in the source and target columns. Apparently the segments table is not being updated with the new identifiers created in the _vertices_pgr table.

An example:

# \d esch.roads
               Table "esch.roads"
 Column |           Type            | Modifiers
--------+---------------------------+-----------
 id     | numeric(10,0)             | not null
 source | smallint                  |
 target | smallint                  |
 geom   | geometry(MultiLineString) |
Indexes:
    "roads_pkey" PRIMARY KEY, btree (id)

# SELECT pgr_createTopology('esch.roads', 0.0001, 'geom', 'id', 'source', 'target', 'true');
NOTICE:  PROCESSING:
NOTICE:  pgr_createTopology('esch.roads',0.0001,'geom','id','source','target','true')
NOTICE:  Performing checks, pelase wait .....
NOTICE:  Creating Topology, Please wait...
NOTICE:  -------------> TOPOLOGY CREATED FOR  8 edges
NOTICE:  Rows with NULL geometry or NULL id: 0
NOTICE:  Vertices table for table esch.roads is: esch.roads_vertices_pgr
NOTICE:  ----------------------------------------------
 pgr_createtopology
--------------------
 OK
(1 row)

# SELECT id FROM esch.roads_vertices_pgr WHERE id NOT IN (SELECT source FROM esch.roads);
 id
----
  2
(1 row)

# SELECT id FROM esch.roads_vertices_pgr WHERE id NOT IN (SELECT target FROM esch.roads);
 id
----
  6
(1 row)

This seems wrong to me, or am I missing something?


texts: Calculating distance from multiple points to one point

texts:
I am attempting to measure distance from the calculated centroid points (using Feature to Point from data management) of every county in the contiguous 48 US states to a single location I have imported with geographic coordinates. I am using the Point Distance tool from the analysis toolbox and I keep getting wacky measurements that do not make any sense. I have used three different county datasets: US census county data, US census cartographic boundary shapefiles county data, and the USA counties layer package from ESRI's arcgis website.

First I attempted to use NAD 1983 contiguous USA Albers equal area projection (at the request of my boss), and realized that my scale for the map was all kinds of wrong (1:300 for the entire contiguous U.S.), and my distance measurements were between 1 and 50 even though they should have been in meters. So I switched to using the North America Equidistant Conic projection, and still was having the same issue. Finally I attempted to use a geographic coordinate system, GCS North American 1983, which fixed the scale issue but I'm still getting the same distance measurements??


texts:
I have some polygon-shapes in qgis-1.7.0 and want to ensure that every polygon is a rectangle. See image below.

Is there any way I can do this in QGIS (or any other f/oss)? Or how can I improve my polygon-drawing to make sure for example that every building consists of rectangles (all angels = 90°)?

Any hints appreciated! Thanks!


texts: I am facing problem when working with Google Earth imagery. When I am digitizing at 1:1000 it just shifts the shapefiles when I even pan the image. Checked for CRS and everything but it doesn't work for me.


texts: Where can I download ArcGIS Service Pack 2?

texts:
I understand SP2 is available for download. http://blogs.esri.com/support/blogs/supportcenter/archive/2011/04/20/arcpad-10-0-service-pack-2-released.aspx

  We are pleased to announce that
  Service Pack 2 for ArcPad 10.0 is now
  available to download at the Esri
  Customer Care website.

Maybe I am missing something but when I log into the ESRI customer care website the Service Pack is not listed.

Anyone?


texts: Export/print QGis: dismatching OSM data on an OSM-layer

texts:
On QGIS (2.10, winXP).
is there any known problem by exporting or printing .osm-data (from JOSM) onto an Openstreetmap-Mapnik-layer?
The imported waypoints always shift from the left less and less to the right-/east side of the map. On the left side, they seem to be good, but on 2.5 Km to the left, they are wrong for about -20 to -100 m - the .osm-data-layer seems to become smaller.
I adjusted the coordinate system for QGis
von QGIS +proj=longlat +datum=WGS84 +no_defs
+proj=longlat +datum=WGS84 +no_defs
the layer of the waypoints is "WGS 84, EPSG 4326"
+proj=longlat +datum=WGS84 +no_defs
the OSM-layer is "WGS 84 / Pseudo Mercator, EPSG 3857"
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
might this mistake caused by mismatching coordinate systems?
But it seems for me, that I can't change the systems?


texts: Changing field types so graduated symbols are available

texts:
I can't seem to make graduated symbols work. When I try to click on my column it doesn't show any options.
I think these is due to those columns being formatted as text and not number. Any ideas how can i work this out?


texts: Raster vertical pixel offset (ALOS AVNIR)

texts:
I am working with ALOS AVNIR-2 level 1a. There is a vertical offset/shift of each and every pixel column. It gives sort of a stripy appearance when zooming in. Does anyone have experience with this and know how to correct it, preferably in envi, arc or gdal.
It does not seem to be software dependent (opened it in ENVI version 4 and 5, ArcGIS, nest, beam)


texts: What is a good plotter that works well with ArcMap?

texts: I will soon be setting up a small office and we will be needing a 36" plotter. I am looking for recommendations for a good plotter that works well with ArcGIS. Speed is not important as long as I can produce good maps for field work and analysis. (Budget = < $5,000 CAD)


texts: What topologic relation means "interior-disjoint"?

texts:
I am using the JTS Topology Suite, and trying to write an expression that means "shape A and shape B are interior disjoint".

I read this documentation page but I don't entirely understand the relations.

The "disjoint" relation (which is the opposite of "intersects") means that the two shapes don't have any point in common, including the boundary.

The "overlaps" relation means that the two shapes have some points in common but no shape is covered by the other. Thus the "not overlaps" relation doesn't take care of the case in which one shape entirely covers the other.

There is a function that defines the "DE-9IM Intersection Matrix", but, I couldn't figure out what matrix implies the "interior disjoint" relation.


texts: ArcGIS Server Object Extenstion REST Endpoint That Supports GET and POST

texts:
I have a Server Object Extension for ArcGIS Server 10.0, that implements IRESTRequestHandler.  The extension contains many individual RestOperations.  By default each RestOperation supports HTTP GETs when using the extension.  In the RestOperation constructor, there is a parameter called postOnly, which is a boolean.  If I set postOnly = true, then the endpoint will only support HTTP POSTs, and will throw an exception if you attempt to use a GET.  I would like to create an endpoint that supports both GET and POST.  If you interrogate some of the supported methods on a standard MapService like Find and Identify, you will notice these support both GET and POST:

http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/identify
http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer/find

Does anyone know of a way to implement this on a custom Server Object Extension?


texts: Nested polygons in QGIS overlap and don't show the real colour gradient

texts:
I am working with nested polygons in QGIS and each of them has got a specific let's say n (numeric) value represented with a given colour gradient.

All my polygons are loaded into one shapefile
If I use Feature blending mode = Normal I miss (don't know why) some of the nested polygons, but colours are ok, see below:

If I use Feature blending mode = Multiply QGIS shows me ALL the polygons BUT for all of them their colours overlap and their areas become darker than the reality, see below:

What I need is to show ALL my polygons with their real (not-overlapped) colour gradient.

Is there a way I can get around this?


texts: QGIS Bookmarks not transfering with files

texts: My students send me QGIS map files that when I open them have been stripped of the bookmarks.  Is there a setting I'm missing or is the normal case.


texts: Georeferencing shapefile using ArcGIS for Desktop?

texts:
I have a shapefile with some polygons that I would like to superimpose to  a raster image, but unfortunately they don't match perfectly due to the reference system the shp has been created with. So, I now want to "georeference" the shp so that all the polygons would match the below raster image. It's not just a matter of moving it, the polygons are a bit deformed compared to the img, so I need to create some control points to tell the program how I would like to move them. I know this can easily be done for rasters, but what about features?

Moreover, I am adding that I know that I can always georeference the raster to the shp, but in this case the raster must stay where it is, I only need to move the shp, as this is supposed to be the second part of the same job, but coming from another source has caused this problem which I have to deal with.

I am using ArcGIS 10.1 for Desktop.


texts: Combining 2+ tables in ModelBuilder into just 1 supertable?

texts:
This might be a stupid question but is it possible to add multiple tables in ArcGIS Modelbuilder together? Just add 2 tables together, like a spreadsheet list without a common factor? I have been googling this for ages and all the methods seems to require some common factor?

For example:

Table 1:
x | 2
z | 5

Table 2:
a | 10
c | 13

Combined Supertable:
x | 2
z | 5
a | 10
c | 13


texts: I brought in a vector shapefile with around a total of 20 different pressures.  I want to make this vector into two different shapefiles.  One for the property "Pressure" of lines that have under 20 pressure and another for the same property that has a pressure of 20 and over.  I know of the Vector > Split Vector Layer command but that splits the one shapefile into 20 different ones.


texts: Reclassify raster based on Z-values in Arcmap 10.4.1

texts: I have a raster which has a Z-value applied to each pixel. I want to reclassify the raster so that any pixels with values below 1 are removed. Is there an ArcMap tool which can perform this task?


texts: OpenLayers freehanddraw forEachFeatureInExtent returns features not in shape

texts:
I am not a newbie.
I am trying to get the features(icons) found inside a closed freehand line.
The approach I am taking is to take the extent returned from the drawing and
use it to get the features(icons) inside the drawing. This works. It includes some extra features(icons).  Is there a better approach to finding featues(icons) inside the freehand drawing that is more accurate? If not I have had some success using the following:

The extent returned with draw.on('drawend',  function (event) {
is used to get features in closed line drawn on map.

vectorSource.forEachFeatureInExtent will have features not inside drawing

var geomA = event.feature.getGeometry();

var bhpoints = geomA.flatCoordinates;

I use these flatCoordinates(xy) to compare to each feature to get a more
accurate decision whether the feature is in or outside of the closed line.

This works but will miss once and a while.

To restate my goal:

I want to find the features(icons) inside of a freehand drawing and the
obvious solution with extents needs better accuracy.


texts:
I imported precinct data from the Harvard Dataverse into Postgres 9.3 using the PostGIS Shapefile Import/Export Manager.

It creates a geom column with a SRID of 0.  I think the data is projected as SRID 4326 as I was able to view a map of it using TileMill when I chose WGS 84.  However, QGis says it is EPSG:4019 so maybe that is the problem.

How can I import it with this SRID?  When I try to change the SRID it fails.

SELECT UpdateGeometrySRID('public', 'precinct', 'geom', 4326);
ERROR:  Geometry has Z dimension but column does not


texts:
We have a problem using the Geoserver version 2.5.2. We are using Geoserver for a line bus fleet monitoring application. We have a data synchronization problem.
I have attached the map that we use to give you an example:

This happens to me if I click on the 1068 bus icon (marked with the red arrow) but the position data has been updated in the view we use (from a server-side service). We do not return any information because the coordinates of that bus have changed. But this creates a serious problem because it seems that there is a malfunction on the map because the bus is not selected and receives no information.

The call we make to the geoserver is the GetFeaturesInfo and we are downloaded png.


texts: Changeing Display cache path permanently in ArcMap?

texts:
I'm looking for a way to change the Cache path of ArcMap permanently.

So from -
   C:\Users\$username$\AppData\Local\ESRI\Local Caches

To a seperate drive -
    E:\Users\$username$\ESRI\Local Caches

I'd imagine a registry entry would control this - I haven't been able to find one.

I'd be happy to run a Python script on startup of ArcMap to change this in a automated way - haven't been able to find any documentation on this.


texts:
Why doesn't the script break if the Version exists?

  Traceback (most recent call last):
    File "C:\Users\jnmiller\Desktop\Create User Version.py", line 30, in
      arcpy.CreateVersion_management(sdeConnection, parentVersion, versionName, "PROTECTED")
    File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\management.py", line 18244, in CreateVersion
      raise e
  ExecuteError: ERROR 001148: Cannot create a version with the specified name.
  Failed to execute (CreateVersion).

# Name: SDE_CreateVersion
# Description: Creates a new SDE version

# Import system modules
import arcpy, time, smtplib, os

# Set local variables
sdeConnection = "Database Connections\NEO_1.sde"
parentVersion = "sde.DEFAULT"
versionName = os.environ.get( "USERNAME" )

fullPath = os.path.join("Windows Credentials: EGAS", versionName)
existingVersion = os.path.join("Version Currently Exists: ",versionName)
createVersion = os.path.join("Version Created: ",versionName)

print("Searching for User Credentials..")
time.sleep(2)
print(fullPath)
time.sleep(3)
print("Checking SDE for existing Credentials..")

for version in arcpy.da.ListVersions(sdeConnection):
    if version.name.split(".")[0] == versionName:
        time.sleep(3)
        print(existingVersion)
        time.sleep(5)
        break
    else:
        # Execute CreateVersion
        arcpy.CreateVersion_management(sdeConnection, parentVersion, versionName, "PROTECTED")
        time.sleep(3)
        print(createVersion)
        time.sleep(5)


texts: Bentley Microstation in comparison to ESRI ArcMap : Filetypes and Capabilities

texts:
I'm curious if there is a side-by-side comparison of Bentley's Microstation v8i (Select series 3) with ESRI's ArcMap 10.1 showing capabilities and filetype acceptance.

I have heard some conversations lately where a Bentley user stated unequivocally that Microstation could do more and open more filetypes natively than ArcMap could.

I'm not a fan of either over the other, but I viewed Microstation as a CAD system with some GIS capabilities whereas ArcMap was a GIS system that could do CAD work. And that compared side-by-side Microstation would be the lower ranking application as ArcMap just had more capabilities. But now I'm wondering if there is something I am unaware of or whether someone else is just quite clueless.


texts: QGIS 'Rasterize vector to raster' failing to create viable raster

texts:
I had a DEM (.asc):

Which was converted to a shapefile point cloud (see attached below) and I'm trying to convert back to raster within QGIS 3.8.2

https://1drv.ms/u/s!AjCedBZJ5Eh4kEKj4AAjcItprRax?e=hPjIX7

The Shapefile loads fine in QGIS, and I can see that points have varying Z values for height

Using Rasterize (vector to raster) I set up with resolution of 641 matching original DEM

But all I get is a black to white colour ramp top to bottom as output:

What am I missing here?


texts:
I'm troubling with save of shapes in db. The shapes are steel sheet cuts, with may be holes, therefore are collections. I decided to solve to my self the problem, to allow the easy exportation between pc. But I cannot achive to detect the mirrored shapes. I wrote some very robust code to detect the equality, but doesn't detect the mirrored copies. The nature of the problem impose to detect it, because if a sheet has a good surface,  these are different shapes. I would wonder if postgis detect mirrored copies as different copies or also treat them as equal?
Ore if somebody can help me with my routine, to detect the equality I wrote an adjacenze matrix of distances between each vertex. If the sum of each row are equal to the new shape the shapes are equal( an extension of a triangle is equal to another if three infos are equal) . But I cannot find a way to define the mirror, if not go throug perverse code and analyze the points.


texts:
I am trying to load a table with over 27mn rows from my spatialite DB using QGIS.
It has seem to hung while loading but I am not sure if it is just taking time to load everything.
Even if it loads, I suspect the performance will be slow.
Are there options to speed up the performance - i.e is there a more efficient format to load large number of datapoints with features? Or should I break my data into smaller chunks? i.e. its 2 years worth of data so instead of combining them i can do each year as a table/layer instead?
Note: The table has a spatial index on it


texts: Passing field name as parameter to script in ModelBuilder?

texts:
I am new to ModelBuilder.

I am running a python script from ModelBuilder that requires as it's parameters a table view and specific field names from that tableview. Essentially the script is looking for information from specific fields without knowing what the user may call the fields.

For example, the script needs
an 'address' field and the user may have identified their address field as 'street' in their table.

I am wanting a dropdown field populated by the tableview that the user can select the correct field from and pass that field name to the script as a parameter. I have tried using
the Get Field Values tool and arcpy.GetParameterAsText(1) returns the value of the field but what I need is the field name itself.

I was thinking there should be a way to create a list and populate it with something like:

fieldList = []
fields = arcpy.ListFields(tablePath)

but I have no idea how I would do this in ModelBuilder and be able to use that list to populate and parameter list the user could then use to select the correct field.


texts: Ogr2Ogr exported this boundary 90 degrees incorrectly

texts:
I've got some GEOGRAPHY data in my Sql Server 2008. I wish to export it as a shapefile (so I can import it into QGIS) and when I do, it's rotated 90 degrees to the left and then flipped 180 degrees to the right.

Command line

C:\>ogr2ogr -f "ESRI Shapefile" "C:\temp\test.shp" "MSSQL:server=dbserver; database=xwing; uid=xwing; pwd=dbpwd"
            -sql "SELECT OriginalBoundary FROM Geography Boundaries where LocationId = 18561"
            -overwrite
            -lco "SHPT=POLYGON"

Boundary is Sql Server 2008

Boundary loaded into QGIS after using Ogr2Ogr

/me sad panda.


texts: How to calculate minimum cell value for a multi band raster?

texts:
I'm a python-qgis beginner.

I use to work with multi-band rasters.
Each raster is made of images coming from satellite sensors acquired in different dates and stacked in order to build a multi temporal dataset.

Is it possible to produce a raster layer in which each pixel has the minimum value that this pixel reaches in all the bands?

I mean, it is possible to produce the "layer of temporal minimum" of the multi-band raster?


texts: Multiple shortest paths with pgRouting using a list as a secondary table

texts:
I have been using a postgis database with pgRouting to solve shortest path problems. I have been using a query like:

select seq, e.geom, pd.cost, pd.id1 as node
from pgr_dijkstra(
'select gid as id, source, target, length as cost
from tl_road',
15489,
12854, false, false)
as pd left join tl_road as e
on pd.id2 = e.gid
order by pd.seq

tl_road is my network. 15489 and 12854 are example start and finish points. This works well.

I have some commuting data of start and finish points and I have merged and found the nearest network node for each point, then imported it into my postgis database. My second table looks like this:

id  RefNo       Origin      Destination
--  --          --          --
1   65277       61780       33025
2   67784       22005       18464
3   65897       257         17407

How can I ammend my pgrouting query to find the shortest path for EVERY ref number on my list using the origin & destination columns?

I have done a lot of googling and have found various solutions but these all seem to be multiple start points to ONE destination or vice versa. I have tried various adaptations of my pgrouting query and usually get the error: "more than one row returned by subquery used as an expression".

Ideally I want to be able to export one result table with one shortest path polyline for each refno id with the refno identifier in the attributes.

At the moment my list is 20 records. I may want to do this for more records in the future but it depends on memory and processing time.


texts: Is Python essential for working in QGIS?

texts: I am a beginner in QGIS. I Generally have a doubt: to work in QGIS is it necessary to know Python? whether i am asking we should  know to write python programs to make certain modules run?


texts: Changing decimal places of X,Y coordinates using proj.4?

texts:
How do I change the accuracy of the output coordinates X,Y in the Proj.4 library (version 4.9)?

In accordance with the manual I tried:

proj +proj=poly -f '%.4f' +R=1

However, for lat=lon=50, the results

'.4f'   '.4f'

indicate that the format has not been recognized (OS: Win10).


texts: Calculating distance between two Linestring using Python

texts:
I have two linestring, each linestring represents a rail road. I am trying to calculate minimum distance between these two linestrings. Which function is efficient to do this?
I have only come across distance function that takes a point and linestring but not 2 linestring.
For e.g.:
l1 = LineString([(42.073407, -87.806245), (42.0752508,-87.8080299)])
l2 = LINESTRING(35.442827 -79.470579, 35.444889 -79.469465, 35.445829 -79.468907, 35.446608 -79.468294, 35.447893 -79.46687)

I need the distance between these two shapely objects.


texts:
I have hundreds of workspace files (FMW) and I would like to import and publish it to FME Server.

Can I do this automatically, not individually by FME Desktop?

If yes, please write me instruction how.


texts: Overwrite a PDF item type with ArcGIS API for Python

texts:
I'm trying to overwrite an existing PDF ArcGIS Online item using ArcGIS API for Python.
props = {
    "type":"PDF",
    "title":file_name,
    "overwrite":True
    }

data = os.path.join(tempFolder, "{}.pdf".format(file_name))

gis.content.add(item_properties=props, data=data)

Returns
RuntimeError: Item already exists


texts: How can I remove a feature's edge without hitting 'Del'

texts:
We use tablet PC with windows 10 and QGIS installed.
We need to modify polygon and polyline features and we need to delete some edges. Normally, we use the 'Del' key of the keyboard but with the virtual keybord of win 10, this key is missing.

We have installed another virtual keyboard but this is not easy to use.
If possible, it would be great getting a button under QGIS  which removed selected edge.

This button exist under the topoedit toolbar but it works only with postgis layer.

Is there a Python script that could do the job?


texts:
I'm a new user of QGIS. I'm georeferencing a scanned image of topographic map. At least 16 GCPs were required for us to locate. However, only 12 can be located and can be identified and be defined based from the coordinates given in the corners, horizontal and vertical borders and tick marks in the middle.

Are there any other ways to obtain define coordinates of other GCPs accurately? Is interpolation possible in this case? If yes, how?


texts: Sub setting shapefile to plot a region and addressing mismatch in ggplot & data

texts:
I am plotting India's district map with data (cases).  First I read data for each district and district shapefile. Later I merge the two on district. Using fortify data I get ggplot. I get a beautiful plot using ggplot.

nr1 <- read.csv("../OUTPUT/working_districts_state.csv")
india.dist <- readShapeSpatial("./shapes2/2011_Dist.shp")
india.dist@data <- merge(india.dist@data,nr1[,c(1,3)],by="DISTRICT")
india.dist@data$id <- rownames(india.dist@data)
ss.fort <- fortify(india.dist)
ss.fort1 <- merge(ss.fort, india.dist@data, by="id")

d <- ggplot() + geom_polygon(data = ss.fort1, aes(x = long, y = lat, group = group, fill=reg.cases)) +coord_equal()

I have two issues:

First, I want to plot a state separately. I sub-set the shapefile and try to plot using similar code. Problem arises when I use fortify command. The command creates an id column with ID extracted from original data i.e. india.dist and not the sub file. So when I merge state.fort with state@data, I do not get proper merge. I am still perplexed as I do not know how to reset ID.

kr <- india.dist[india.dist$ST_NM== "Kerala",]
state <- kr
state@data$id <- rownames(state@data)
state.fort <- fortify(state)
state.fort1 <- merge(state.fort, state@data, by="id")
unique(state.fort[,c(4:7)])
     hole piece  id group
1    FALSE     1   7   7.1
363  FALSE     1 173 173.1
719  FALSE     1 232 232.1
1384 FALSE     1 277 277.1
1960  TRUE     2 277 277.2
2013 FALSE     1 288 288.1
2561 FALSE     1 314 314.1
3018 FALSE     1 320 320.1
3349 FALSE     1 321 321.1
3753 FALSE     1 358 358.1
4144 FALSE     1 422 422.1
4735 FALSE     1 436 436.1
5204 FALSE     1 567 567.1
5663 FALSE     1 571 571.1
6292 FALSE     1 610 610.1

state@data
             DISTRICT  ST_NM ST_CEN_CD DT_CEN_CD censuscode reg.cases id
1           Alappuzha Kerala        32        11        598       259  1
2           Ernakulam Kerala        32         8        595       769  2
3              Idukki Kerala        32         9        596        42  3
4              Kannur Kerala        32         2        589        50  4
5           Kasaragod Kerala        32         1        588        38  5
6              Kollam Kerala        32        13        600       381  6
7            Kottayam Kerala        32        10        597        63  7
8           Kozhikode Kerala        32         4        591       237  8
9          Malappuram Kerala        32         5        592        95  9
10           Palakkad Kerala        32         6        593        31 10
11     Pathanamthitta Kerala        32        12        599        30 11
12 Thiruvananthapuram Kerala        32        14        601        97 12
13           Thrissur Kerala        32         7        594        99 13
14            Wayanad Kerala        32         3        590        39 14

Secondly, I suspect that the ggplot is not matching data (cases) to district properly. It is mismatching the data. I do not know how to check that such mismatch is not happening and if it is happening then how to fix it. I think this issue also involves handling ID in shapefile.


texts: Units of measurement for vector grid in QGIS?

texts:
Can anybody tell me as to what the units of measurement are (e.g. meter, kilometer, miles etc) in vector grid menu of QGIS (see attached image)?

I want to make cells of certain sizes in meter or kilometer (e.g. 50 meter, 100 meter, 500 meters etc.) but my coordinates are geographic, which is making it difficult for me to specify the cell size in meters or kilometers.


texts: An algorithm that could find a random cycle in a graph

texts: I am looking for an algorithm that could find a random cycle in a graph from a node while that cycle is traversing around another nodes (area). For example, from the green star on the left of the image, finds a random cycle that goes around the red-star nodes.


texts: Auto batch for all of the layers in my mxd

texts:
I have a mxd file with more than 800 layers, and I want to convert all of these layers to KML files. I have tried to use the batch option in the "Layer to KML" tool, but I have to set there every layer manually, which is crazy.
Is there any way to auto-fill all of the layers in my mxd file into the batch table? There is any other way to convert multiple layers into a KML files? (One KML file for each layer).


texts: netCDF raster layers not loading in QGIS

texts:
I have worked with netCDF files in QGIS before, but since I installed QGIS 3.6 they won't load ("myfile.nc" is not an admitted raster data source), nor do they load in QGIS 2.18 anymore. I find this particularly strange. I am working under Windows 10. It appears both versions of QGIS are using GDAL 2.2.4. I don't know if that is the problem but when I run gdalinfo --formats on the OSGeo4W Shell, netCDF does show up:

netCDF -raster,vector- (rw+s): Network Common Data Format

I was looking for a way to downgrade GDAL in case this was the problem but did not find anything feasible at my level. Any help would be welcome.


texts: Convert .ply and .lin files to any file format readable by QGIS?

texts:
My issue:
I'm trying to open .ply and .lin files in QGIS but they do not appear to be supported, (although this might also have to do with the fact that there are actual issues with my data - see below).

My question:
How can I convert the .ply and .lin files to something/anything that QGIS can read?

Here are my troubleshooting steps so far:

First of all, here are a few descriptions of the file types:

PLY File Format - Description 1
http://paulbourke.net/dataformats/ply/

PLY File Format - Description 2
https://en.wikipedia.org/wiki/PLY_(file_format)

PLY files can be in ASCII or BINARY form. I believe mine are in ASCII form, according to the README text document that came with the datasets.

One major place I'm running into trouble is that:
Apparently PLY and LIN files require a header, but my data does not appear to contain headers, per both me not seeing any headers when viewing the data in Atom and an "unable to locate headers" error I received when trying to open the files in MeshLab, (open source software that allows viewing and converting of PLY files).

Similar SE Question:
I've also seen this question regarding converting .ply to .las, but the solution seems a bit convoluted for my purposes, and without any headers I don't think anything is going to be able to read my files at this point anyway.

Next steps / new question:
When I view the files in Atom, I can see all of the attributes. I can actually add the headers listed in the README text doc, adding them to the PLY and LIN files in Atom, and then I can save them out as CSVs, which I can then open and view in Excel. What I see in the PLY file is all of the attributes of each feature. What I see in the LIN file is all of the geometry (lat/longs). Both files do have unique IDs that can be used to join up both of these files.

However, since these are polygons, and not points, it's not an easy enough join with each record having only 1 pair of lat/longs. Some records have a handful of lat/long pairs, and some have over a dozen. In many cases one feature has multiple rows of lat/long pairs, which I'm interpreting as multi-part features.

So I have all of these data, but I can't figure out a way to manually convert these CSVs into some kind of spatial format. The problem is I have all of the data I need, I just don't know how to make it spatial.

Snapshot of PLY csv:

Snapshot of LIN csv:


texts:
I tried to modify spikeRemover from Schmidt & Krüger. My goal it's to obtain, for each polygon, the position of the spikes. I don't need to remove it.
For some reason, that I can't figure out, my script fails and I don't obtain any position.
Probably there is a problem at the union of different points but I don't know how to manage it in another way.

    CREATE OR REPLACE FUNCTION ST_SpikeFinder (geometry, angle double precision)
     returns geometry as
     $body$
     DECLARE
     ingeom alias for $1;
     angle  alias for $2;
     lineusp geometry;
     newgeom geometry;
     numpoints integer;
     point_id integer;
     result_pnt geometry;

    begin
            -- input geometry or rather set as default for the output
        newgeom := ingeom;

        -- check polygon
        if (select st_geometrytype(ingeom)) = 'ST_Polygon' then
            if (select st_numinteriorrings(ingeom)) = 0 then

            lineusp := st_boundary(ingeom) as line;
                    numpoints := st_numpoints(lineusp);
            point_id := 0;

                -- the geometry passes pointwisely
            while (point_id <= numpoints) loop
                        result_pnt := ST_Collect(result_pnt, st_pointn(lineusp, point_id));
                    -- the check of the angle at the current point of a spike including the special case, that it is the first point.
                    if (select abs(pi() - abs(st_azimuth(st_pointn(lineusp, case when point_id= 1 then st_numpoints(lineusp) - 1 else point_id - 1 end),
                                st_pointn(lineusp, point_id)) - st_azimuth(st_pointn(lineusp, point_id), st_pointn(lineusp, point_id + 1))))) <= angle then

        (probably the problem is here)      --->        result_pnt := ST_Union(result_pnt, st_pointn(lineusp, point_id));

end if;
            point_id = point_id + 1;
            end loop;
        end if;
        end if;
            return result_pnt;
     end;
     $body$
       language 'plpgsql' volatile;


texts: How do I turn off some formats in GeoServer's WMS?

texts:
I have published a couple of layers on my Geoserver.

To be complaint with the T&C of my data provider, I need to block the access of some of those layers as Vector data. Hence I have already blocked the WFS services on these Layers.

I came to know that even on the WMS service, you can get Vector Data, if you use GeoRSS and some of the other formats.

How do I block GeoRSS, SVG & UTFGrid as formats in the WMS Service?


texts: Check if error due to no spatial reference of extracted raster?

texts: I have raster file without spatial reference and I extracted my area of interest by shapefile that has spatial reference. Finally, I got extracted raster with spatial reference. Can I assume the raster taken the spatial reference from the shapefile during extraction process?


texts: How to calculate field based on another fields

texts:
I´m trying to calculate field based on another fields. I´ve got this attribute table.
Example: If field "Smisenost" is S than FRST in new field.
Example 2: If field "Smisenost" is NULL, if field "KULTURAKOD" is NULL, if field "ZpusobyVyu" is NULL, if "DruhPozemku" is 7 than RNGE.
Field "DruhPozemku" has always value, but another ones not. So I need script for calculate value from values that you can see in image.

It should work by this key.


texts: How can i delete multiple nodes, for a polygon, with Quantum GIS?

texts:
I know how to select a single node and delete it.
I can also click-drag to select MULTIPLE nodes. But if i click delete, only one node (from the higlighted list) is removed.
Can I delete all the highlighted nodes?
UPDATE:
Here's a quick YouTube clip, showing how it's not working. U can see that when I select 2 or 3 nodes, only ONE is deleted.


texts: Error whilsted using QGIS 3.4.4 r.series batch processing

texts:
I'm having some trouble using the Batch Processing to perform r.series, Could anyone explain what I'm doing wrong? I'm using QGIS 3.4.4.

Output Log:

  Processing algorithm 1/1… Algorithm r.series starting… Input
  parameters: {'-n': False, 'GRASS_RASTER_FORMAT_META': 'XML',
  'GRASS_RASTER_FORMAT_OPT': 'GTiff', 'GRASS_REGION_CELLSIZE_PARAMETER':
  10.0, 'GRASS_REGION_PARAMETER': None, 'input': 'D:\LocalDataStore\Pre-Processing\DecloudedRasters\18NVM\Band_02\Declouded_RT_T18NVM_20151221T153112_B02.tif;D:\LocalDataStore\Pre-Processing\DecloudedRasters\18NVM\Band_02\Declouded_RT_T18NVM_20151224T153612_B02.tif',
  'method': [2], 'output': , 'quantile': 'float',
  'range': '-1,1e+08', 'weights': '1'}

  g.region n=0.0 s=0.0 e=0.0 w=0.0 res=10.0 r.series method="median"
  quantile="float" weights="1" range="-1,1e+08"
  output=outputf220a810b6eb44f884953602e7378fcf --overwrite g.region
  raster=outputf220a810b6eb44f884953602e7378fcf r.out.gdal -t -m
  input="outputf220a810b6eb44f884953602e7378fcf" output="D:\TEST_3.tif"
  format="GTiff" createopt="GTiff" metaopt="XML" --overwrite Cleaning up
  temporary files... access: No such file or directory Traceback (most
  recent call last): File
  "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py",
  line 413, in processAlgorithm Grass7Utils.executeGrass(self.commands,
  feedback, self.outputCommands) File
  "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 386, in executeGrass feedback.reportError(line.strip())
  TypeError: reportError(self, error: str, fatalError: bool): not enough
  arguments

  Batch execution completed in 1.62 seconds


texts: Can Cartodb map show only United States?

texts: I have made my map with cartodb, showing where certain work is going on throughout the United States. I want to send it out, but would prefer the map to only show the United States--whereas now it shows Canada and Mexico as well. I am using the free cartodb.


texts: Adding PostGIS layer from QGIS Python plugin

texts:
I try to develop a plugin to add postgis layers to qgis from a plugin as part of learning python.

The code for adding a postgis layer is taken from the cookbook and work fine if I run it from the python console within qgis (If button OK is pressed in dialog box result == 1 and should add my layer to qgis).

But if I run it from my plugin it throws the message NameError: name 'QgsDataSourceURI' is not defined. Why do I get that error running it from the plugin?

Is there a difference how I add a layer from a plugin/inside a function vs from the python console?

 def run(self):
        """Run method that performs all the real work"""

        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec_()
        # See if OK was pressed
        if result == 1:
            # Do something useful here - delete the line containing pass and
            # substitute with your code.
            uri = QgsDataSourceURI()
            uri.setConnection("localhost", "5432", "test", "postgres", "postgres")
            #set database schema, table name, geometry column and optionaly subset(WHERE clause)
            uri.setDataSource ("basic", "cities", "geom")
            # Defining the layer name and layer type for QGIS?
            vlayer=QgsVectorLayer (uri .uri() ,"cities","postgres")