texts: Why AddJoin creates 'one to many' style of output?

texts:
What I am trying to do is to join feature class (2 records/polygons for a test purpose) with excel table(70 records but only 10 would have the same join value, 5 to each record/polygon) using join field created and calculated with use of data from other columns.

What I am getting is a feature class with 10 records (each record created 5 times)

What I want to have is a feature with 2 records + fields from the joining table.

What I am doing excel spreadsheet (with 3000+ records) converted to table (gdb) -->from table i am exporting only 70 with the same 'Option_Name'--> adding joining field --> calculating field (it is text, 150char eg: 15FGS00012/1)

then:

original feature class converted to work feature class (gdb) - adding joining field --> calculating field (same format as in table)

So far so good, I am getting join fields properly calculated, but when I am exporting (or appending, or printing) joining ouput I am always getting that 10 polygons.

What I am doing wrong?

import arcpy

Report = r"P:\GIS\FGS Datasets\Data\application_16102016_copy.xls\application__contr$"

output_gdb = r"P:\GIS\Datasets\Data\Data.gdb"
Report_table = r"Report_table"
optionTable = r"S_TABLE"
expression = '"Option_Name" = \'Rural Woods\''

field2 = "GIS_Join"
field2_type1 = "TEXT"
field2_concat1 = "Application_Submitted_Reference"

SDR = r"J:\ArcSDE Connections\GIS.sde"
FC_raw = "GIS.S_RURAL"
FC_FGS = "S_RURAL"
FC_concat1 = "SCHEME_REF"
FC_concat2 = "FGS_OPTION"
FC_featureLayer = "S_RURAL_layer"
expression2 = r"Rural Woods"

final_gdb = r"P:\GIS\Datasets\Data\Data_Final.gdb"
final_feature = final_gdb + "\\" + FC_FGS

arcpy.env.workspace = output_gdb
arcpy.env.overwriteOutput = True
workspace = arcpy.env.workspace

# export Report.xls to output_gdb as a table
arcpy.TableToTable_conversion(Report, output_gdb, Report_table)
print Report_table, " created"

arcpy.TableSelect_analysis(Report_table, optionTable, expression)
print optionTable, " created"

arcpy.AddField_management(Report_table, field2, field2_type1, "", "")

cur1 = arcpy.UpdateCursor(Report_table)
for row in cur1:
    if row.Option_Name == "Rural Woods":
        field2_concat3 = "1"
        expression3 = row.getValue(field2_concat1) + r"/" + field2_concat3
        row.setValue(field2, expression3)
        cur1.updateRow(row)
print field2, " in ", Report_table, " calculated"

arcpy.FeatureClassToFeatureClass_conversion(SDR + "\\" + FC_raw, output_gdb, FC_FGS)
print FC_FGS, " created"
arcpy.AddField_management(FC_FGS, field2, field2_type1, "", "")

cur2 = arcpy.UpdateCursor(FC_FGS)
for row2 in cur2:
    if row2.FEATDESC == " Publ Access - Rural":
        field1 = "1"
        expression5 = row2.getValue(FC_concat1) + r"/" + field1
        row2.setValue(field2, expression5)
        cur2.updateRow(row2)
print field2, " in ", FC_FGS, " calculated"

arcpy.MakeFeatureLayer_management(FC_FGS, FC_featureLayer)

arcpy.AddJoin_management(FC_featureLayer, field2, Report_table, field2, "KEEP_ALL")

del workspace
arcpy.env.workspace = final_gdb
arcpy.env.overwriteOutput = True
workspace = arcpy.env.workspace

arcpy.DeleteRows_management(final_feature)
print "All data (rows) deleted from the previous version of ", final_feature

arcpy.Append_management(FC_featureLayer, final_feature)

del Report, cur1, cur2


texts: Getting attribute table into Excel spreadsheet?

texts:
How can I get an attribute table into an excel spreadsheet?

I don't have the tool for that because I am using ArcGIS Desktop 10.0.


texts: Converting Historical Railroad Systems Data from KML?

texts:
I'm writing a paper for an undergraduate economic history class and I'm interested in North American railways in the 1800's. I've found historical GIS data for the period here. The railroad data are KML files.

My idea is to analyze the railways using network/graph theory. Each station would be a node/vertex, and the edges would be the rail lines between each station. I could then use measures like closeness centrality to analyze the network structure.

To do so I would need to convert the KML files to basically a spreadsheet (i.e .csv) format giving the distance between each connected station, with row x column y giving the length of the rail line between station x and station x (an adjacency matrix).

Is it feasible to perform such a conversion?

My undergraduate degree is in economics and mathematics but I have experience with Python from data analysis and computer science courses. I don't have any direct GIS experience.


texts: What is the search mechanism in "join attributes by location"?

texts:
While running join attributes by location, we can either choose one to many or one to one mapping. In One-to-one, it assigns the first located features value.
 What is meant by first located feature? How is this locating or searching carried out?


texts: how can I specify a popup to show only for a specific layer in openLayers 3?

texts:
I am using the function to define popup contents in a map, and right now, it createa a popup for every feature I click on the map, no matter what layer it is in. I would like this popup to appear only for features within one layer. I'm sure there is a way I can specify this in the function, but I'm not sure how. I tried using this for line 6 if (feature && layer === layerVector) { but it gave me an error saying that layer is undefined. layerVector is the variable name for the layer I want to use the popup for. Can anyone help? Here is the code:

map.on('click', function(evt) {
  var feature = map.forEachFeatureAtPixel(evt.pixel,
      function(feature, layer) {
        return feature;
      });
  if (feature) {
    popup.setPosition(evt.coordinate);
      content.innerHTML = '<b>Parcel No: </b>'+
      feature.get('parcel')+'<br>'+
      '<b>Owner: </b>'+
      feature.get('primowner')
  }
});


texts: Using Boolean Parameter in QGIS Graphic Modeler

texts:
I am trying to create an if/else condition with a boolean parameter in QGIS graphic modeler. The parameter will be used inside Field Calculator2 algorithm.By default the parameter is False.

I am using the boolean parameter to create an if/else condition inside the formula of the Field Calculator2 Algorithm. See below. The model is working, but the condition always give result for False output even when I make Boolean Parameter True.

Is there something that I am doing wrong in the if/else formula here?

if( 'drainedsoilsleaveuncheckedifnotsure' = True,replace("GDCodeTemp", '/D', ''),replace("GDCodeTemp", map('A/', '', 'B/', '', 'C/', '')))


texts: How to use Fields in the Field Calculator?

texts:
I'm using the Field Calculator in ArcMap and I need to create a unique ID for every storm drain in my county.
An ID Should look something like this: 16-I-003

The first number is the municipal number which is in the column/field titled "Munic"

The letter is using the letter in the column/field titled "Point"

The last number is simply just 1 to however many drains there are in a municipality.

So far I have:

rec=0
def autoIncrement()
pStart=1
pInterval=1
 if(rec==0):
  rec=pStart
 else:
  rec=rec+pInterval
return "16-I-" '{0:03}'.format(rec)

So you can see that I have manually been typing in the municipal number, the letter, and the hyphens. But I would like to pull the first number and letter from the fields: Munic and Point so I don't have to manually type them in each time it changes.  I was hoping I could just put !Munic! and !Point! in my return statement before the .format(rec) or somehow work it in there.

I'm a beginner when it comes to python and ArcMap, so please dumb things down a little. The simpler the better!


texts: Introduce year instead of whole date in date type field using ArcGIS Desktop?

texts:
I have a field with numbers in it (concentration of chemical compounds that were presented in aquatic environments) and I'd like to combine them with a field which would have date type capabilities but as a unique year (e.g. 2013, 2008 etc.) instead of a whole date (2/5/2013, 5/8/2008 etc.).

The year data represent the years of chemical analysis of water samples.

The general concept of this work is to provide me finally the information of the concentration and the year of analysis of a certain compound.

Also I'd like to have the ability to select concentrations by year.

For example: before the year of 2004 what was the concentrations of these substances?

How can I achieve the above?


texts: How can I connect my point pattern shapefile with my polygon window in R?

texts:
I would like to do some spatial statistics in R. I have a point layer called '1983' that represents the distribution of one species, and a polygon layer called 'SA-R' that represents the study area.

I followed the steps from the pdf "Handling shapefiles in the spatstat package":

  library(spatstat)

spatstat 1.36-0       (nickname: ‘Intense Scrutiny’)
For an introduction to spatstat, type ‘beginner’

  library(sp)

  library(maptools)

Checking rgeos availability: FALSE
        Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
        which has a restricted licence. It is disabled by default;
        to enable gpclib, type gpclibPermit()

  Year1 <- readShapeSpatial("1983.shp")

  Studyarea <- readShapeSpatial("SA-R.shp")

Field name: ‘opp mÂ²’ changed to: ‘opp.mÂ²’
Field name: ‘opp ha’ changed to: ‘opp.ha’

  as.owin(Studyarea)

window: polygonal boundary
enclosing rectangle: [226675.38, 234961.33] x [21233.66, 32693.92] units

  as.ppp(Year1)

marked planar point pattern: 7 points
Mark variables: ID, Year
window: rectangle = [227669.54, 229542.88] x [22341.009, 30989.048] units

  Data1 <- Year1[Studyarea]

Error in x[j] : invalid subscript type 'S4'

I can plot each of my shapefiles in R, but do not manage to plot my point layer within my study area polygon. Does someone know what's wrong with my last command?
I have already checked, and both of my shapefile layers are in the same projection system.

Your syntax works until this step:

  w -> as.owin(Studyarea)
  Error in as.owin.default(Studyarea) : Can't interpret W as a window

For some reason, it did work like this:

  library(spatstat)
  library(maptools)
  library(sp)
  Year1 <- readShapeSpatial("1983.shp")
  as.ppp(Year1)

marked planar point pattern: 7 points
Mark variables: ID, Year
window: rectangle = [227669.54, 229542.88] x [22341.009, 30989.048] units

  Studyarea <- readShapeSpatial("SA-R.shp")

Field name: ‘opp mÂ²’ changed to: ‘opp.mÂ²’
Field name: ‘opp ha’ changed to: ‘opp.ha’

  w <- as.owin(Studyarea)

  Year1 <- ppp(coordinates(Year1)[,1],coordinates(Year1)[,2], window=w, marks=Year1@data$MyVar)

  plot(Year1)


texts: Running bounding-box queries with Google Maps

texts:
I am working on a google map application, to display the items that are listed in the selected rectangular area. Below image give you some understanding.

Based on the co-ordinates (topLeft,bottomRight), i am running bounding box query to find the items in  mongo db. But it always returns the this error when running the box queries

 Error: $box needs 2 fields (bottomLeft,topRight)

I know that i have passed topLeft,bottomRight points instead the mongo expected one (bottomLeft,topRight).

How do I calculate the bottomLeft,topRight values based on the co-ordinates from Google Maps?

This is the actual value i got trouble with

bound_box = [[80.22050443383796, 13.051049734963444], [80.2393871853028, 13.043106313697518]]]

and ruby code for finding bounding box values in mongo db

items = Item.where(:loc => {"$within" => {"$box" => [bound_box]}},
                              :category => BSON::ObjectId(cat), :sub_category => BSON::ObjectId(sub)).order_by(:created_at => :desc).skip(@current_page * 5).limit(5).paginate

this is the line giving problem

items = Item.where(:loc => {"$within" => {"$box" => [bound_box]}})

which is similar to the mongodb query

 db.items.find({"loc" : {"$within" : {"$box" : bound_box}}})


texts: Shapefile to SVG with ArcMap (while retaining country attributes)

texts:
Does anybody know how to export the Natural Earth 1:110m shapefile to SVG with ArcMap while retaining the country attributes?

From SVG I'm gonna try to convert it again over to RaphaelJS to achieve something like:

http://backspace.com/mapapp/javascript_world/

http://backspace.com/mapapp/javascript_world/js/world_570.js

But I'm new to all of this so I'm not sure how to go about it. Maybe there's no way but to manually add the country codes for the JS?


texts: Can't get OpenLayers 3 to display data from QGIS Server WMS service

texts:
I have set up a QGIS server and I can get to the data it spits out via a WMS connection in QGIS desktop.  Now I'm trying to get the same service to display in an OpenLayers 3 map.

Here's the JavaScript I'm trying to use:

var map;

function init(){
    map = new ol.Map({
        target:'map',
            renderer:'canvas',
            view: new ol.View({
            projection: 'EPSG:2966',
            center:[2809414,1031652],
            zoom:5
        })
    });

    var newLayer = new ol.layer.Image({
        source: new ol.source.ImageWMS({
            url: 'http://localhost/qgis/qgis_mapserv.fcgi.exe',
            params: {'LAYERS': 'Lateral', 'FORMAT':'image/png',
           'VERSION': '1.3.1'}
        })
   });

    map.addLayer(newLayer);

Nothing is displayed on the page, and I'm getting a error in the JavaScript console:

TypeError: e is undefined ol.js:56

If I need to post the html let me know, but there isn't much to it.


texts: How to configure mapserver correctly on cherokee?

texts:
I am attempting to configure mapserver using fastcgi on cherokee and here's what I have got at the moment:-

Unfortunately, that doesn't seem to be enough as I somehow need to have mapserv running as a process.  But I don't see such an option for mapserver in its docs.

Any suggestion how I may configure mapserver correctly on cherokee?


texts: How to query multiple features from Geoserver WFS?

texts:
Are there any methods to query several features by id other than one entry by a single url request?

The scenario is I store around 200 point layers in PostgreGIS. Each maintains 100k entries GPS points. My boss wants to let them "move" in the browser. Obviously it's impossible to query all entries at once. So I want to query like 1k entries for one layer in every query. So basically I want to query the features by their id, like "starID<=ID<=endID". But it seems like by WFS url query, it can only query one entry every time by the featureID=ID. I'm also thinking about using the ordinary Postgre other then PostgreGIS, but just in case there might be spatial queries in the future.


texts: Measuring Geodesic Distance in ArcGIS for Desktop?

texts:
I'm working with only a basic/student license for ArcGIS 10.1 on an ocean mapping project.

I've created a grid that covers all ocean area, capturing various characteristics within that grid (the grid is now in point data).

I now need to add the distance from each point to the edge of various polygons (shipping lanes, countries) but need geodesic distance because of the large distance covered.

I read that if you keep the data in a geographic coordinate system, unprojected, the distances generated through various measuring tools in ArcGIS automatically measure in geodesic.

Is this true, and if not, is there a way to force geodesic measurements in various tools?


texts: Web Service level IP restriction with ESRI Web Adaptors?

texts:
In order to restrict service usage when sharing a REST endpoint I'm looking for a way to control access of service based on an IP range with the ESRI web adapter.  I know that this can be done with a reverse proxy but they only documentation I have found for doing this with ESRI's web adapter is cross domain restriction.

Not all of my services require the same amount of security, some are public while others are not, so the solution needs to target individual services with different security needs.

If this can't be accomplished with the above method, what is the best method for restricting service REST endpoint access. Using the ESRI server OOB users and roles restrictions does not prevent username and password sharing, when sharing url's.  I have to be certain unauthorized users are unable to access specific published services.

Server is AWS based on Linux (RHEL)


texts: Is there any statistical measure/tool that compares two different sets of location configurations in R or ArcGIS/QGIS?

texts:
Here is what I'd intend to do;

I have a given set of facility locations as shown below;

And based on a predefined one optimization criterion, I obtained the following configuration (and there are many other which resulted in different location configurations);

I understand that whatever the optimization criterion I used, I can simply look at the criterion value for both and conclude that the optimized network is "better" (lower for minimization problem and higher for maximization problem). But I am in need of another statistical measure that will tell how similar the new network (i.e., 2nd figure) is when compared to the existing network (i.e., 1st figure). Below is the figure showing both configurations;

There are several different optimized configurations obtained through utilizing different optimization criteria and the intent here is to quantitatively assess which optimization criterion would provide a solution (configuration that is most similar to the existing network.

I am thinking that calculating an average of shortest distance of all pairs could do the job but this may require some kind of a matching algorithm that I am not aware of.

Can anyone here suggest an existing package or tool that can do the job described above?


texts: Creating\rendering a heatmap not showing up

texts:
I am having a problem rendering a heatmap. I am sure the heatmap is created but it doesnt display when I replace or add a renderer.

  var TF = new FeatureLayer({
    url:"https://services3.arcgis.com/1pxU2hJU9ZszJDcX/ArcGIS/rest/services/Young_Adult_Treatment_Facilities/FeatureServer/0",
  })

  //apply HM to layer
  function apply(result){
    console.log("applying")
   TF.renderer = result.renderer;
   };

  var heatmapParams = {
    layer: TF,
    basemap: map.basemap
  };

document.getElementById("butt").addEventListener("click", heatmapRendererCreator.createRenderer(heatmapParams).then(apply));

Whenever I call the apply function to create a renderer the point layer disappears.

Here's a codepen:

https://codepen.io/segev-salman/pen/eaEbRv


texts:
Using ArcGIS 10 I want to create a fishnet with its bottom left coordinate = (-450000,-1401771.01064457)
I then want 71 rows above from this point and 79 columns right.
Based on 36000m by 36000m cells.

I am not sure if the one coordinate i have is the origin or the y a-axis coordinate.

How do I find the other coordinate so I can create the fishnet?


texts: Show overlay on top when zoom out in OpenLayers 4.6.5?

texts:
I have a tile layer, three vector layers.  The tile layer has lots of overlays in it. I want to put vector layers on top when zooming out the map using OpenLayers 4.6.5. I tried to apply setLayerIndex and raiseLayer, but it says it is not a function of OlMap. Is there a way to get this done?

Layer instantiations:

initializeMap() {
    this.source = new OlVector({});
        this.clusterSource = new OlCluster({
          distance: this.distance,
          source: this.source
        })

        this.layer = new OlVectorLayer({
          style: (f, r) => this.styleFunction(f, r),
          source: this.clusterSource,
        });

        this.dockLayer = new OlVectorLayer({
          style: (f, r) => this.styleDockFunction(f, r),
          source: new OlVector({}),
        });

        this.countLayer = new OlVectorLayer({
          style: (f, r) => this.styleLabelFunction(f, r),
          source: new OlVector({})
        });

        let tileLayer = new OlTileLayer({
          source: new OlXYZ({
            crossOrigin: null,
            url: 'https://cartodb-basemaps-{a-c}.global.ssl.fastly.net/light_all/{z}/{x}/{y}.png'
          })
        });

        let view = new OlView({
          center: OlProj.fromLonLat([5.73271, 58.96564]),
          zoom: mapZoomLevel,
        });

        this.map = new OlMap({
          interactions: OlInteraction.defaults().extend([
            new PinchZoom()
          ]),
          layers: [tileLayer, this.layer, this.dockLayer, this.countLayer],
          view: view
        });
}

Marker additions to map follows:

createMarker(pos, txt, color1, txt2, color2, dockingStationId) {
    let elem = this.marker.nativeElement.cloneNode(true);
    elem.querySelector("#marker-txt").innerHTML = txt;
    elem.querySelector("#docking-station").innerHTML = dockingStationId;
    elem.querySelector("#marker-txt").style.backgroundColor = color1;
    elem.querySelector("#marker-txt2").innerHTML = txt2;
    elem.querySelector("#marker-txt2").style.backgroundColor = color2;

    var marker = new OlOverlay({
      position: pos,
      positioning: 'center-center',
      element: elem
    });
    this.map.addOverlay(marker);

  }


texts: Printing highly detailed and/or complex maps on large format paper?

texts:
I don't know if this is the right place for this question but, since it involves maps created with ArcMap, hopefully someone met this same problem.
I want to print a highly detailed map of the whole town of Florence (Italy) comprising of a topographic basemap and about a hundred thousands points representing delivery points for postal service (see picture).

I'm going to use an HP Designjet T520 36in plotter and want to print it on an A0 format paper. I've tried both printing it directly from ArcGis (though I'm using a remote version of it) and exporting it as PDF and print it from Adobe Reader (PDF size is around 30 Mb). After I made sure all page and printing settings were ok, I pressed the print bottom. Apparently the printer doesn't like it since I just see a printing bar lingering for some seconds on 0% and for an imperceptible istant on 100%, but nothing happens. The bar just disappears. I suppose it depends on the size of my map, but I'm not sure, and anyway I don't know how I could go around that. I've searched a lot through the web but hasn't found anything helpful yet. So, how would you print detailed and complex maps? Have you ever met a problem similar to mine?


texts: Trying to figure out ArcGIS API's `mapScale` variable and how it maps to OpenStreetMap's zoom level

texts:
I'm using ArcGIS' /export endpoint to get an .svg of a few map layers given a specific bounding box. However, when I apply the .svg as an image overlay to OpenStreetMap's world map using Leaflet, the overlay is scaled slightly too small.

I figured out that the mapScale parameter to the /export endpoint could be changed to alter the scale. However, there's no way for me to know what "scale value" my Leaflet map uses. As far as I can tell, Leaflet only knows about zoom levels and not "scale values".

Are there any resources for me to figure out how to find the proper scale value when exporting from ArcGIS?

Here's how I'm adding the svg layer:

var imglayer = 'layer.svg';
var imgbounds = [
    L.CRS.EPSG3857.unproject(L.point([-10791685.4014,4999593.1461])),
    L.CRS.EPSG3857.unproject(L.point([-8301672.7680,7822259.7266]))
];
L.imageOverlay(imglayer, imgbounds).addTo(mymap);


texts: Picking coefficient in Inverse Distance Weighted method of QGIS?

texts:
I'm using QGIS 1.8.0 to interpolate weather data. This is, I have a weather station located in say, A and C. With the Inverted Distance Weighting or Triangular Interpolation methods, I can know/interpolate information in B. When using the IDW method, I have to choose a coefficient.

Does anyone know in which range the coefficients should be?

Say from 1 to 10 or .1 to .99 or 100 to 1000?

Any ideas?


texts:
I have a Raster, which covers only the southern part of the country, like this:

I wish to have a Raster, which covers the entire country (including the North part)? The northern part should have the nodata values.

How do I do this using GDAL or QGIS?

( I do not have ArcGIS, so I can't use the answer given here: Increasing raster extent by filling extent with no data using ArcGIS Desktop without Spatial Analyst?)

I tried using gdalwarp's -te parameter; It was successful in reducing the extent, but did not extend the extent to the Northern part.


texts: Visualizing overlay of 6 point data sets using ArcMap

texts:
I have a visualization question. Using ArcMap 10.6, I would like to overlay six point datasets in a way that synergies and differences can be quickly seen by people who have no idea about GIS applications and are not use to work with digital maps.
My current idea would be to simply use different symbols and colors for the six datasets and display this on a map. This representation could be combined with a simple density analysis (kernel density) (see figure). I am not sure if correlations would be understandable.

Nonetheless, I would like to ask if there are more "understandable" ways to display six point datasets overlaid so that the differences and similarities can be quickly identified by anyone?
Information about the points: The points represent potential locations for new trees that were freely chosen in a group discussion (six group discussion results). Each point record could contain between 3 to 50 points. First, all groups should be compared with each other.


texts:
To reproject a table into different CRS is a very common practice in geospatial analysis. However, I am still struggling to reproject my table in PostGIS, after reading loads of Q&A related to this topic.
Here I got a table with a set of points. The original SRID is 28356
SELECT Find_SRID('public', 'pts', 'geom');

I'd like to reproject the table to 4326 for further analysis.
Firstly, I used ST_Transform, but failed with this error message:
ERROR:  Geometry SRID (4326) does not match column SRID (28356)
UPDATE pts SET geom = ST_Transform(geom,4326);

Then I tried ALTER table but got the same error: ERROR:  Geometry SRID (4326) does not match column SRID (28356)
ALTER TABLE pts
  ALTER COLUMN geom
  TYPE Geometry(Point, 4326)

Finally, the UpdateGeometrySRID function seems running correctly and the Find_SRID query return 4326 as expected.
However, when loading the modified layer to the map canvas, they are not in the correct location.
SELECT UpdateGeometrySRID('public', 'pts', 'geom', 4326);
SELECT Find_SRID('public', 'pts', 'geom');

Is there anything wrong with my query?


texts: contingency tables or cross tabulate in Earth Engine

texts:
Working on Google Earth Engine and still finding this challenging. I have two images, the first one is an image of landcover categories and the second one land management units. I want to cross these categorical images to obtain basically a count of pixels of each landcover in each unit.

I guess this is called a cross-tabulation, a contingency table or a tabulate operation (in ArcGIS).  I have noted the tutorial on the reducer to get zonal statistics but this is not quite what I need https://developers.google.com/earth-engine/reducers_grouping.

Is this possible to implement in Earth Engine? I would be surprised if not!

Vectorizing one image to do the analysis seems like a very lousy approach.


texts: Creating Adjacency Matrix using SpatialLinesDataFrame and SpatialPointsDataFrame in R?

texts:
I am trying to create an adjacency matrix of power plants that are connected via power lines. The plant data is point data (~8000 observations) with longitude and latitude measurements. The line data must have been mapped/drawn since it does not include specific location measurements.

I created 2 feature classes (one for power plants and another for power lines) and created a geometric network in ArcGIS in which many of the points snapped to the lines (there are ~7000 points that now intersect with lines). I believe the snapped/modified data was saved to the feature classes.

I imported the 2 feature classes into R using readOGR (as part of the rgdal library) and they converted into SpatialLinesDataFrame (for the power lines) and SpatialPointsDataFrame (for the power plants) classes. The import worked correctly since the plot appears to be exactly like that in ArcGIS.

The next step is to build an adjacency matrix (~7000x7000) where a 1 indicates a connection between two powerplants (data in the SpatialPointsDataFrame) via the intersecting powerlines (data in the SpatialLinesDataFrame) and a 0 otherwise. There are around 1000 plants that do not intersect with lines so these will either not be included in the matrix or will not have any connections.

Packages, such as shp2graph, can create adjacency matrices, but appear to only do so with SpatialLinesDataFrame (adjacency matrix that indicates where the lines intersect). I am wondering if there is a package/method that can make this adjacency matrix of points that are connected through lines. It is assumed that two plants are still connected even if they are through a series of intersecting lines.


texts: How does a GeoTIFF file store elevation data?

texts:
I'm trying to understand how elevation data is stored in GeoTIFF files, when used for Digital Elevation/Surface Models. My very basic understanding is that in a GeoTIFF file (or at least those produced by Fusion's dtm2tif function), each pixel has a value associated with it for each band, and that these values are 0 if no data is present, and between 1 and 255 for data.
Is there a piece of metadata in the GeoTIFF that indicates how to scale between 1-255 to the elevation data in the elevation model? Is there an easy way of accessing this scaling in python?
Cheers
Ben


texts: Choosing Equal Area Projection and Standard Parallels for Indonesia?

texts:
I am trying to do area computation stuff over whole of Indonesia and as such, am using the following projection:

http://spatialreference.org/ref/esri/asia-south-albers-equal-area-conic/

To make this projection work over Indonesia, I am following the guidelines give at the following link:

https://geonet.esri.com/thread/30524

The recommended modifications change the following:

central meridian: 115.0 (from 125)

standard parallel 1: 2.0 (from 7)

standard parallel 2: -7.0 (from -32)

latitude of origin: 0.0 (from -15)

I think central meridian and latitude of origin do not impact distortions in any manner. However, why is it that the standard parallels should be 2 and -7? The entirety of Indonesia does not seem to be lying within these latitudes. From my inspection of the map, the entirety of Indonesia would lie somewhere between 6 and -11. Wouldn't these be better choices? How does a change of standard parallels/secants impact the computations?


texts: Transform equirectangular simple TIFF to Mercator projection with gdal_translate / gdalwarp

texts:
I'm trying to transform this shaded map of Venus, a 2000x1000 TIFF in what I believe is equirectangular projection, without GIS metadata, to Mercator.

I first tried

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3395 \
  venus_shaded.tif venus_shaded-mercator.tif

but this fails with

ERROR 1: The transformation is already "north up" or a transformation
between pixel/line and georeferenced coordinates cannot be computed for
venus_shaded_cyl.tif. There is no affine transformation and no GCPs.
Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this
check.

Simply adding SRC_METHOD=NO_GEOTRANSFORM is clearly not the answer, as it produces garbage:

I then tried using gdal_translate to add GCPs:

gdal_translate \
  -gcp 0 0 0 90 \
  -gcp 0 333 0 30 \
  -gcp 0 666 0 -30 \
  -gcp 0 999 0 -90 \
  -gcp 666 0 120 90 \
  -gcp 666 333 120 30 \
  -gcp 666 666 120 -30 \
  -gcp 666 999 120 -90 \
  -gcp 1333 0 -120 90 \
  -gcp 1333 333 -120 30 \
  -gcp 1333 666 -120 -30 \
  -gcp 1333 999 -120 -90 \
  venus_shaded.tif venus_shaded-gcp.tif

(and then running the same gdalwarp command on the new file), but the results make it clear I don't know what I'm doing.

I don't know whether the problem is that I'm misunderstanding the format of the -gdp option to gdal_translate, misunderstanding how to apply it to pixel coordinates, misunderstanding the whole concept of GCPs, or all of the above. It's also possible EPSG:4326 is the wrong source coordinate system.

Regardless, I feel like this should be a pretty straighforward operation, and there's probably a better way to do it. What is that way?


texts:
I have two Layer in my Project

1) Voronoi Layer which is constructed on point Layer

2) Building Layers which is represented as Polygon.

The goal is to find the nearest node for every polygon (building).

I am using QGIS 3.0 and using python console within QGIS to automate this.

I am highly confused in below two approaches: NN join and Nearest Neighbour.

I am thinking to find the centroid of each polygon and then find the nearest node using the Distance Matrix Tool.

Small Green dots are Nodes and Polygon are Buildings.


texts: Why can't I merge DEM rasters in QGIS 3.10?

texts: I recently downloaded QGIS 3.10 and am trying to merge several raster DEM .tif files via the Ratser tab > miscellaneous > merge. I tried several output data types like Float32, Int16, and Byte, and high and low compression profiles but they all yield the same error message: "The following layers were not correctly generated (insert layer name). You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm." I am not sure where to find this information or what it means. Do I need to install a plugin to run it? What can I do to make this command run successfully?


texts: How to sort a field according to value range?

texts:
I have to check an attribute table field (MY ID NO) of polygon feature against a tiff image of parcel map where fields are written as 12-58,65,87,2565-89999,58648.

Example: Parcel map tiff(image) contains plot no (MY ID NO in poly feature's attribute table) against each plot, each tiff has a corner box that describe the plot numbers included in this tiff map. Now I have to check is all the opposite plot no (collected from tiff image file, and also mentioned in the corner-box describing included plot no) is populated in the MY ID NO field.

Now I want  a field calculator script (I think it is easy and better way to check, suggestions also preferred) to populate "Sorted Range" field. I need only "Sorted Range" field. My table as below:


texts:
I'm working on a project that involves using linear interpolation to fit lines to
lat/lng points. I'm pretty new at scientific programming related stuff, and I was wondering if people could steer me in the right direction.What is a good tutorial for linear interpolation in Python, preferably using GIS data?


texts: We use Saex extension in ArcGIS for parcel updating. In some computers we get the error Error Occurred in Processing. I have attached the screenshot for reference.


texts: Why does gdalwarp return incorrect corner coordinates for my geotiff?

texts:
I have a png source image that's in Lambert Azimuthal Equal Area and I have the long/lat extent. So I ran gdal_translate to create tie points and convert to gtif. Then I ran gdalwarp to reproject to Web Mercator. gdalinfo shows that the bounding coordinates are ok for the gdal_translate output, but not for the gdalwarp output. Que pasa?

Here's my gdal_translate operation on the source file:

gdal_translate -of GTiff -a_srs "+proj=laea +lat_0=40. +lon_0=-105 +x_0=0 +y_0=0 +a=6370997 +rf=298.257" ^
-gcp 0 0 -119.344674 47.148697 ^
-gcp 525 0 -90.655326 47.148697 ^
-gcp 525 431 -93.609897 31.438244 ^
-gcp 0 431 -116.390103 31.438244 ^
-a_ullr -119.344674 47.148697 -93.609897 31.438244 ^
src.png tmp.gtif

(I included the bounding box long/lats for good measure but it seems to make no difference in the outcome of gdalwarp.)

Here's my gdalwarp operation on the gdal_translate output (tmp.gtif):

gdalwarp -of GTiff -s_srs "+proj=laea +lat_0=40. +lon_0=-105 +x_0=0 +y_0=0 +a=6370997 +rf=298.257" -t_srs EPSG:3857 tmp.gtif final.gtif

Here's the metadata from 'gdalinfo tmp.gtif':

Driver: GTiff/GeoTIFF
Files: tmp.gtif
Size is 525, 431
Coordinate System is `'
Origin = (-119.344674000000000,47.148696999999999)
Pixel Size = (0.049018622857143,-0.036451167053364)
Image Structure Metadata:
  INTERLEAVE=PIXEL
**Corner Coordinates:
Upper Left  (-119.3446740,  47.1486970)
Lower Left  (-119.3446740,  31.4382440)
Upper Right ( -93.6098970,  47.1486970)
Lower Right ( -93.6098970,  31.4382440)
Center      (-106.4772855,  39.2934705)**
Band 1 Block=525x3 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=525x3 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=525x3 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=525x3 Type=Byte, ColorInterp=Alpha

Corner coordinates look good. But here's the metadata from 'gdalinfo final.gtif':

Driver: GTiff/GeoTIFF
Files: final.gtif
Size is 579, 355
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+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"],
    AUTHORITY["EPSG","3857"]]
Origin = (-11688702.286414048000000,4866004.054033085700000)
Pixel Size = (0.057992644981843,-0.057992644981843)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
**Corner Coordinates:
Upper Left  (-11688702.286, 4866004.054) (105d 0' 5.04"W, 40d 0' 1.53"N)
Lower Left  (-11688702.286, 4865983.467) (105d 0' 5.04"W, 40d 0' 1.02"N)
Upper Right (-11688668.709, 4866004.054) (105d 0' 3.95"W, 40d 0' 1.53"N)
Lower Right (-11688668.709, 4865983.467) (105d 0' 3.95"W, 40d 0' 1.02"N)
Center      (-11688685.498, 4865993.760) (105d 0' 4.49"W, 40d 0' 1.28"N)**
Band 1 Block=579x3 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=579x3 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=579x3 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=579x3 Type=Byte, ColorInterp=Alpha

Also notice that the pixel size is suspect in the final image - both dx and dy are the same.


texts: Cutting a raster file utilising a vector polyline

texts:
I’m on 3.16 QGIS GRASS... It hasn’t got gdalwarp.
Simple issue, I have a raster file which includes two regions: land and sea. I have a vector polyline that marks the division between land and sea. I want to keep the sea element of the raster, and cut, mask, bounded out, loose or otherwise not see the land element.
Seems I have to turn the polyline into a polygon first? (Or drop to MS paint, manually trace, re geo reference, and reload with sellotape and tippex ..)?


texts: Need to clip raster based on field name in ArcGIS 9.3

texts:
Is there a way around for me to batch clip a raster using a field name to create output a raster layer ? Much like using Split tool but it is for a vector. How about that for the raster side ?

Going to do that for slope and elevation by 24k quad. I am using ArcGIS 9.3.1 and I am aware that the 10 has a raster split.


texts: Release License from ArcGIS License Administrator

texts: I want to release two machines using my Network Analyst extension from license server administrator. The problem is that those two machines no longer exists on our domain and it is impossible to release license from that machine itself. The machines are listed on Borrowed tab on the license server administrator but not on checkedOut tab. What can I do about this issue?


texts:
I am beginner and new for app RTKGPS+ in smartphone. I am interest to use this program to use my smart phone as the rover by use GPS signal from itself. And take the data from base station to this program for my survey onsite with accurate required less than 1m.

Can I do this with this program and how I can do that?

Can I also get information about manual for this program?


texts: UTF-8 CSV import ok but after processing with utf-8 shapefile accent replaced by '?'

texts:
I m using QGIS 2.18

I export data from SQL Developer to a csv file with ";" as separator and with UTF-8 encoding.

Then, I have imported the csv with UTF-8 encoding. Accents were imported right.

But after processing a join by attributes between CSV and a Shapefile with .cpg "UTF-8" the accents disappear:

What's the matter? Is it from ogr2ogr sub "processing" that does not support the accents?


texts:
I would like to convert .shp into .gpx.

Using GIS, I open my shapefile and I tried to "save as" gpx but it gives an OGR error:

  creation of field ObjectId failed (OGR error: Field of name 'ObjectId'
  is not supported in GPX schema. Use GPX_USE_EXTENSIONS creation option
  to allow use of the  element).

How can I fix this error?


texts:
I'm currently trying to pull neighborhood data from Zillow's Neighborhood page (http://www.zillow.com/howto/api/neighborhood-boundaries.htm) and then load it onto a leaflet map, openstreetmap.  I created a "Neighborhood" table and added the data into PostGIS using shp2pgsql.  Everything works great, I can query the data and get back GeoJSON using ST_AsGeoJSON, however the query is really slow.  I'm considering how I could speed up this process as it looks like ST_AsGeoJSON is the primary culprit.  This table will not be updated that often, so I'm wondering if it would make sense to simply store the GeoJSON in another column within Neighborhood and then return that instead of converting the geom column every page refresh.

Here's an example of my Neighborhood table.

CREATE TABLE "neighborhood" (gid serial,
"state" varchar(2),
"county" varchar(43),
"city" varchar(64),
"name" varchar(64),
"regionid" numeric);
ALTER TABLE "neighborhood" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','neighborhood','geom','0','MULTIPOLYGON',2);


texts: Seeking comprehensive information on how to create beautiful and artistic cartography?

texts:
In many circumstances I experienced how powerful maps can be because a key decision maker just enjoys looking at a beautiful representation of facts instead of reading 10 pages of text. I have a background in geography, and we  learned GIS software at the university, however the courses never discussed in depth aspects of beauty and aesthetics.

I would like to learn, however, more about what makes a map really beautiful and I'm searching for literature/blogs/websites/courses etc. Aspects like color choice, fonts, layouting for different formats (web/print) should be discussed. Basically it goes into the direction of graphic design but in the field of GIS and based on the possibilites of GIS software (preferentially QGIS).

I forgot to tell that (as a start) I already saw Kenneth Field and his wonderfull blog. Kenneth also recently published a comprehensive book that is based on ArcGIS.


texts: Field Calculator - conditional function and string variable?

texts:
I know this question has been asked before in different fashions, but I can't seem to get my python script to work despite looking at all the relevant questions I could find.

I am simply trying to populate an attribute field based on the conditions of another attribute field.  Specifically I have one field that labels various projects under different headings (example: easement, BPS, donated etc.)

I want to calculate the values in a new field to group these into one of three headings (Donated Easement, Purchased Easement, BPS).  However I can't get the python to work and I would rather not do it by hand.  Here is the code I ran just to try to get one label type recalculated.

pre-logic:

def reLabel(typeLabel):
  if typeLabel == 'BPS Easement':
    return 'BPS'
  else:
    return 'Other'

Label2:

reLabel(!Type_Label!)


texts: Label Style is converting decimal String names to Number

texts:
The feature.attributes.fid returns a String value like "9.3210". But on the map, this value is being converted to some number type and on map is showing the value 9.321 (without the 0).

var contextResource = {
    getLabel: function (feature)
    {
        return feature.attributes.fid;
    }
};

var resourceStyle = new OpenLayers.Style({
     label: '${getLabel}'
},{ context: contextResource });

I've tried to make things like:

return JSON.stringfy(feature.attributes.fid); On the map it shows with quotes, like "9.3210"
return String(feature.attributes.fid); This still shows 9.321 (without the 0)

Any suggestions to display 9.3210 (without quotes, and with 0) on map?


texts: Splitting census areas based on another polygon shapefile

texts:
I have one shapefile contains census areas with population and density information and the other shapefile is just an administrative zone file. There may be multiple census areas within each zone and the reverse may true as well.

I would like to calculate the density based on the zone information and but the density information is in the census area. If possible, when the census area is split by the zone boundary, I would like the density to be proportional to the area. For example, if 75% of the census area goes into the Zone A, the density for the census area should be divided so that it has 75% of the density.

I've attached a picture of the two shapefiles.

I run ArcGIS 10.3 on a basic license, so the split tool won’t be an option for me.

Is this doable?


texts:
I have two different shapefiles:

point shapefile with many elements
area shapefile with many regions

I need to create a rule where I need to colour differently the points in some areas, but the areas are in another shapefile

I use QGIS 2.6.1


texts:
GeoServer will throw an error (parserError Error: getJson was not called) for cca 8 out of 14 AJAX requests if I use AJAX in for loop. I have to call AJAX with 500ms interval - then it works without errors. My problematic AJAX:

for (var i = 1; i <= 14; i++) {
      $.ajax({
        url: url,
        dataType: 'jsonp',
        jsonpCallback: 'getJson',
        success: function(result) {
            console.log(result)
        }

    });
}

When I use Python I get all responses right:

import requests
i = 1
while i <= 14:
    r = requests.get("http://localhost:8080/geoserver/d99/ows?service=WFS&version=2.0.0&request=GetFeature&typeName=d99:day_time_" + str(i) + "&sortBy=time&outputFormat=text/javascript&format_options=callback:getJson")
    i += 1
    print(r.text)

I run GeoServer on my notebook (localhost) and GeoServer is not installed as service (I run it manually). Does GeoServer have any limit for WFS requests per second?


texts: How to generate a tiled map from a spherical mercator SVG?

texts:
I'm not a GIS person; this is really just a side task, so I probably have terminology/strategy incorrect; I would ask that you please bear with me. I have the following requirements for a Qt-based C++ application that I'm working on:

Display a tiled map (e.g. OpenStreetMaps). I have this working well.
Allow the overlay of other vector data on top of the map (as a test, I have a shapefile of world boundaries I got from Natural Earth.

Preprocessing of the vector data is possible in my application, so I decided to convert the EPSG:4326-projected Natural Earth shapefile to an SVG file (which I can render most easily in my application). I have done so using Kartograph. Its configuration looks like this:

{
   "layers": [{
       "src": "ne_10m_admin_0_countries.shp",
       "class": "world_boundaries"
   }],
   "proj": {
       "id": "proj4",
       "projstr": "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs "
   },
   "bounds": {
       "mode": "bbox",
       "data": [-180, -85.0511, 180, 85.0511]
   }
}

I believe this should take the input shapefile, reproject it to Spherical Mercator (which should match the projection of the OSM tiles), and generate a corresponding SVG. This all "works," but the projection of the SVG file seems to be a bit off, it seems. When I overlay the SVG features atop the OSM image tiles, they don't quite line up, like this:

The black one-pixel-wide lines are rendered from the SVG. I'm wondering if there's something conceptually that I'm doing wrong in the preprocessing steps. I think that I'm rendering the SVG tiles correctly (I'm treating the SVG document as a square image that is 360 degrees wide and 170.1022 degrees tall and just slicing out the square portion of that image that corresponds to each tile that I need). However, since GIS isn't my job, I'm sure my ignorance could be leading me astray somewhere.


texts: GeoserverTerrainProvider plugin displays faults in elevation tiles

texts:
I try to display elevation data using the triple Geoserver/Openlayer3/Cesium. It seems I have the same problem than https://stackoverflow.com/questions/25936968 and for now, I can not get a working solution!

I have done nearly the same gdal processing to obtain my picture pyramid and whatever the GeoserverTerrainProvider solution I use (BIL or with SLD) I get sharp cut elevation at tile edges (see: http://postimg.org/image/wy0j1ijnp/)

In the geoserver layer preview, my elevation image as bad tile border in 'Tiling: Tiled' mode but is ok in 'Single tile' mode. This may be the cause of my problems: ie. GeoserverTerrainProvider is badly requesting tiles from geoserver?

Or did I fail at generating the pyramid geotiff?


texts:
I am relatively new to this whole geotopic. I have to deal with the following scenario:

From a userperspective, I do three steps

I draw a straight line across a riverbed
I get a resultset from the server back
A diagram is drawn where the x-axis has relative distanceinformation from where the points are along the line and the y-axis gives me altitude information of the riverbed

As a programmer, I get GeoJSON which contains Multiline String information (because the data is segmented?).

The overall structure is simply an array of two dimensional arrays

[
 [[long, lat, alt], ...$segment1],
 [$segment2],
 ...
]

To retrieve the distance information (where in relation to a starting point is a given point?) I pick the first segment's first point triple as a referenec point.

From there I calculate successive the distances for the following points.

As a result, I would expect something like this:

segment1: [0, y1], [1, y2] ... [n, yn]
segment2: [n+1, yn+1] ...

What I really get is something like:

segment1: [0, y1], [1, y2] ... [n, yn]
segment10: [2, y10] ...

Which I find quite confusing.

I assumed the segments to be ordered somehow, so that I get a rising monotonously result.


texts:
Is it possible to prevent certain results based on "class" and/or "types"?  For example, I do not want any results that have the following class or types:

class: natural, highway
type: peak, footway, tertiary

Currently the results for this query includes MANY "motorways".  I only want places (houses, cities) to be returned:

q=1505+North+S&addressdetails=1&exclude_place_ids=206651397%2C232991588%2C69754855&countrycodes=us&viewbox=-105.5274544234%2C39.774601254671%2C-105.0274544234%2C40.274601254671&accept-language=en-US%2Cen%3Bq%3D0.9

My other options include: 1)removing the data from the Nominatim database. 2)  removing results from the resultset after the query has executed.  I don't love either of these other options - but willing to experiment!


texts: GPS post processing - Differential Correction

texts:
I have a large set of collected GPS points that I'm pulling into Arcgis from SQL and creating an xy event layer to then map them. The fields for the points are as follows: recorded, lat, long, altitude, speed, hAccuracy, vAccuracy.

Overall the points are fairly accurate, but of course there are a few cases where the accuracy is off by quite a bit. I have worked with differential correction before, but it was with the same software/hardware (Trimble/Pathfinder) used to collect the data.

Is it possible to correct the points after I have created a shapefile? If not, is it possible to correct with the Pathfinder software them before I bring them into GIS? I also realize that if in fact I can correct them I will have to pull use multiple base files as the GPS points collected are over multiple dates.

Any ideas or help is greatly appreciated!


texts: Change font size of label with data defined expression in QGIS?

texts:
I am trying to change the font size of the feature labels based on attributes from a field.  I can't find an example of how to do this.

I believe I can do this by setting up a data defined expression in the Labels>Text>Size expression but I am not sure how to formulate the actual expression.

I am looking for an example of how this expression would be set up.


texts: Are ArcGIS' UTM zone designations simplified?

texts:
I was recently given a document showing the UTM zones as they appear over the whole Earth. Each zone is 6 degrees longitude by 8 degrees latitude. A numbering system proceeds from 1 to 60 from -180 degrees longitude eastwards, and a lettering system proceeds northward from "C" at -80 degrees latitude. All of this is familiar to me.

In ArcGIS, UTM zones are given with the same numbering system based on longitude, but appear to have done away with the latitudinal lettering system, instead using just "N" or "S", meaning North or South of the Equator. Was this simply done to reduce the number of options users would have to trawl through, or is there a more theory-driven line of reasoning here? Are we losing accuracy by doing away with those zone subdivisions?


texts: Seeking ArcGIS focal statistics alternative in open source Python

texts:
I am struggling to do the pixel aggregating of raster in open-source python as similar to the ArcGIS Focal statistics function does, I would like to make a 5 x 5 rectangular window on which the program's function will calculate the mean of the center pixel using neighbor pixels falling inside the defined window. My input raster values are in float format 0 - 1. Please can anyone suggest, a possible way to do it in python?
I tried the below code, it's not working
import time
import glob
import os
import gdal
import osr
import numpy as np

start_time_script = time.clock()

path_ras=r'D:\Firm_SM\F1A/'

for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')):
    rasterfile_name=str(rasterfile[rasterfile.find('IMG'):rasterfile.find('.tif')])

print ('Processing:'+ ' ' + str(rasterfile_name))

ds = gdal.Open(rasterfile,gdal.GA_ReadOnly)
ds_xform = ds.GetGeoTransform()

print (ds_xform)

ds_driver = gdal.GetDriverByName('Gtiff')
srs = osr.SpatialReference()
#srs.ImportFromEPSG(4726)

ds_array = ds.ReadAsArray()

sz = ds_array.itemsize

print ('This is the size of the neighbourhood:' + ' ' + str(sz))

h,w = ds_array.shape

print ('This is the size of the Array:' + ' ' + str(h) + ' ' + str(w))

bh, bw = 5,5

shape = (h/bh, w/bw, bh, bw)

print ('This is the new shape of the Array:' + ' ' + str(shape))

strides = sz*np.array([w*bh,bw,w,1])

blocks = np.lib.stride_tricks.as_strided(ds_array,shape=shape,strides=strides)

resized_array = ds_driver.Create(rasterfile_name + '_resized_to_52m.tif',shape[1],shape[0],1,gdal.GDT_Float32)
resized_array.SetGeoTransform((ds_xform[0],ds_xform[1]*2,ds_xform[2],ds_xform[3],ds_xform[4],ds_xform[5]*2))
resized_array.SetProjection(srs.ExportToWkt())
band = resized_array.GetRasterBand(1)

zero_array = np.zeros([shape[0],shape[1]],dtype=np.float32)

print ('I start calculations using neighbourhood')
start_time_blocks = time.clock()

for i in xrange(len(blocks)):
    for j in xrange(len(blocks[i])):

        zero_array[i][j] = np.mean(blocks[i][j])

print ('I finished calculations and I am going to write the new array')

band.WriteArray(zero_array)

end_time_blocks = time.clock() - start_time_blocks

print ('Image Processed for:' + ' ' + str(end_time_blocks) + 'seconds' + '\n')

end_time = time.clock() - start_time_script
print ('Program ran for: ' + str(end_time) + 'seconds')

MOdified code based on @Neprin suggestion, however, I would like to modify it based on my file structure,
Please help on this
import numpy as np
import gdal
import cv2
import matplotlib.pyplot as plt
import seaborn as sns

img = gdal.Open('20180305.tif').ReadAsArray() # i have multiple raster i.e.20180305, 20180306, 20180305 so on
 # i want put give the path of folder where i kept my input raster
img2 = np.zeros(np.array(img.shape) + 10)
img2[5:-5,5:-5] = img  # fix edge interpolation
kernel = np.ones((5,5),np.float32)
dst = cv2.filter2D(img2,-1,kernel)/25

# Save the output raster in same name as input with projection


texts: How to create layer in GeoServer and not publish it?

texts:
I'm currently using GeoServer REST API in my Django/Python application to make such routine operations as creating workspaces, raster stores and layers. This is how it looks like in my code:

# create "my_workspace" workspace
requests.post("http://localhost:8085/geoserver/rest/workspaces",
              data="<workspace><name>my_workspace</name></workspace>",
              headers={"Content-Type": "text/xml"},
              auth=("admin", "password"))

# create raster store
requests.post("http://localhost:8085/geoserver/rest/workspaces/my_workspace/coveragestores?configure=all",
data="<coverageStore><name>my_layer_store</name><workspace>my_workspace</workspace><enabled>true</enabled><type>ImagePyramid</type><url>file:C:/path_to_data/"</url></coverageStore>",
headers={"Content-Type": "text/xml"}, auth=("admin", "password"))

# and finally, let's create a layer:
requests.post("http://localhost:8085/geoserver/rest/workspaces/my_workspace/coveragestores/my_layer_store/coverages?recalculate=nativebbox",
data="<coverage><name>my_layer</name><title>my_layer</title><srs>EPSG:4326</srs></coverage>",
headers={"Content-Type": "text/xml"}, auth=("admin", "password"))

In the GeoServer UI I see a separate operation to Publish a layer. However, in my case, when I create a layer I see, that I can watch how it looks like, say, in OpenLayers. So, it seems like my layer gets published automatically. And I do not like it. I want to have a separate REST API operation to publish a layer, and I do not want my layer to be automatically published as a result of layer creation operation. So, how can I achieve this?


texts:
I am trying to import ArcPy module into spyder and it does not work.

  import arcpy Traceback (most recent call last):

  File "", line 1, in
      import arcpy

  ModuleNotFoundError: No module named 'arcpy'

I am using Anaconda 64 bit, Python 3.7, ArcGIS 10.3.1.

I have tried the following, nothing works:

Downloaded Anaconda 32 bit, Python 2.7. I am not sure if have changed everything correctly though. Do I have to change anything else after I installed Anaconda 32-bit with Python 2.7?
Change spyder Python interpreter in presettings to C:\Python27\ArcGIS10.3\python.exe. When I restart Spyder the console is broken.
Add ";%PYTHON_HOME%" to my environment variables in Windows advanced system settings

What else should I try?

I did not see ArcPy importing error after setting path? before.

My problem is solved now.

The following helped:

Installing Python Anaconda 32-bit with python 2.7 and activating python 2.7 in the Anaconda Prompt (https://docs.anaconda.com/anaconda/user-guide/tasks/switch-environment/)
Step 2&3 on the following page:
https://support.esri.com/en/technical-article/000012334
Adding this to my code:
sys.path.append("C:\Python27\ArcGIS10.3\Lib\site-packages")
sys.path.append(r'C:\Program Files (x86)\ArcGIS\Desktop10.3\arcpy')
sys.path.append(r'C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcToolbox\Scripts')
sys.path.append(r'C:\Program Files (x86)\ArcGIS\Desktop10.3\bin')
sys.path.append("C:\Python27\ArcGIS10.3\Lib")


texts:
I want to display the content of a local TMS (a tile map, stored in a directory structure Z/X/Y.png) in QGIS.

This is for quick check and comparison.
I am not willing to set WMTS webserver for this purpose.

I tried to load this kind of file in QGIS but I got only a white or stripped map

  <GDAL_WMS>
  <Service name="TMS">
    <ServerUrl>file:///vdata/myTMS/${z}/${x}/${y}.png</ServerUrl>
     <Format>png</Format>
     <Transparent>TRUE</Transparent>
  </Service>
  <DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>
    <TileLevel>11</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>bottom</YOrigin>
  </DataWindow>
  <Projection>EPSG:3857</Projection>
  <BlockSizeX>256</BlockSizeX>
  <BlockSizeY>256</BlockSizeY>
   <BandsCount>4</BandsCount>
 </GDAL_WMS>

I am running QGIS 2.14


texts: Joining a table with an If/Then Statement

texts:
I want to join two tables in ArcMap, but not every feature has the same common attribute in the two tables I want to join. For example, "Feature 1" is in both Table 1 and Table 2 and both tables have the same address for Feature 1, so that would be a good field to base the join on. However, "Feature 2" is in both Table 1 and Table 2, but only Table 1 has the address for Feature 2. Both tables have the same Company name for Feature 2, though, so Company Name would be the ideal field to base the join on.

Is there a way to join two tables with 1 field, if available, if not, then join using a field 2?


texts: Question about ArcGIS vs. QGIS for kernel density estimations

texts: It appears from the QGIS documentation that the Heatmap Plugin can generate a kernel density surface from a point map and make use of the attributes of the original points to influence this surface via the "Use weight from field" option. Does anyone know whether this plugin in QGIS performs the same function as the "Kernel Density" tool in the Spatial Analyst toolbox in ArcGIS?  (in the ArcGIS tool, the "population field" can be used as a weighting option--I'm assuming that if this weight is zero, the resulting surface is essentially a "heat map"?).


texts: Adding new records to attachment match table from select number of features from feature class using ArcPy?

texts:
I am seeking some ArcPy/Python help. I use the mobile app GIS Pro to collect data in the field, which allows me to export shapefiles containing a field called Photos, containing a list of names of photos associated with each feature. The list appears like this:

I am trying to create a tool that will

Import the shapefile into a project geodatabase
Add the features from the shapefile to a "master feature class"
Update an attachment match table to include the photos associated with the new features
Add the new attachments

I am struggling with step 3. There needs to be a new entry in the match table for each photograph. The match table contains the fields REL_OBJECTID & ATT_NAME. The REL_OBJECTID needs to coincide with the OID assigned to the feature after it has been merged into the "master feature class".


texts: Understanding EOL while scanning string literal from Reclassify in ArcPy?


texts: How can I create a Voronoi diagram on gmplot?

texts:
So I have some coordinates that I'm plotting on GoogleMapPlotter.
Now I've been trying to create a Voronoi diagram for some of my coordinates but I haven't been able to find a library that supports this.
There is geovoronoi but that is mainly used with GeoPandas.
Is there a way to create a Voronoi diagram on gmplot?


texts: Printing portion of a map for multiples locations in QGIS?

texts:
I need to generate .png or .pdf files with a portion of the map that Im currently working on for each object in a shape file. Im working on QGIS3.

I have a shape file with multiple points wich represent the objects that Im interested in. I need to generate one map for each point in the .shp. The map has to be centered on the point and must have an specific scale (1:500). It also has to include other layers that are part of the project.

I was thinking that this could be done with de PyQGIS console because the shape file has more than 500 points and using the Layout Manager would take infinite time.

The image below represent what I need. The point of interest is C143. The map is centered on that point and it also shows objects from other layers that are close to C143 (like text and C145).

Thanks!


texts: Access traffic data using Google Maps API

texts:
The following screenshot shows a traffic filter, adopted in the new Google Maps.

Can we implement the similar function using Google Maps API?

It looks like that google.maps.TrafficLayer can only display the live traffic data.


texts: Unsplitting lines with ArcGIS tool 'Unsplit Lines'

texts:
I have tried to use the ArcGIS tool 'Unsplit lines' to get one polyline from multiple touching lines. Somehow, it does not work with some features. Have a look at the pic nr 1.
 In the pic you see one door from a 'door layer'. With this (and other doors) the unsplitting did only partially work, as it is still separated into four lines. However, I had a look at the vertices, and they share relevant vertices. Still, they don't merge. When I had a look at the vertices the connection is only start node with start node, or end node with end node. As such, I concluded, that with the Unsplit tool only start node and end node will form a connection? I wish I could include an adequate picture showing the vertices of all lines. However, in ArcGIS you can view the vertices of only one feature at a time.
Now my question: Is there another possibility to get the door to become one single feature? Dissolve is not an option as I have many many doors in the layer and, unfortunately, the line features do not have a field saying to which door they belong. As such I cannot dissolve per field. The result is, therefore, a layer in which all doors belond to one single multipart feature, which is not what I need. I need each door being a single feature, this should be done by using some kind of automation as there are many many doors.

Of course I had more than one look in the internet to find answers and tried these tools in ArcGIS: Unsplit lines, Dissolve (not an option), Merge (as this needs selection, this is also not an option). I tried aggregate, but it throws an error, besides that it would give some problems as some doors overlap.

Accordingly I would probably need a tool, which merges line features which share a node without that the sharin nodes need to have the properties start node-end node. Do some of you know how to do this?


texts: Resampling to 0.5° resolution in Earth Engine

texts:
I'm trying to reproject a 500m MODIS image into 0.5° longlat for comparison with climate data. I've followed the instructions for resampling here, but I'm confused about reducing the resolution. I suspect that I'm supposed to change the scale value in reproject(crs, crsTransform, scale), but this value is supposed to be set in meters, which are not consistent across degrees. How can I best convert to a lower resolution in degrees rather than in meters?

What I have so far is:

var modis = ee.ImageCollection("MODIS/006/MCD12Q1").select('LC_Type1');
var region = ee.Geometry.Rectangle(-127.18, 19.39, -62.75, 51.29);

var modisyr = ee.Image(modis.first()).select('LC_Type1');
var modisProjection = modisyr.projection();
var modisResolution = modisyr.projection().nominalScale();
print('MODIS projection:', modisProjection);
print('MODIS resolution:', modisResolution);

// Get the forest cover data at MODIS scale and projection.
var modisReproj = modisyr
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 1024
    })
    // Request the data at the scale and projection of the MODIS image.
    .reproject({
      crs: 'EPSG:4326'
    });

var proj2 = modisReproj.projection();
print('MODIS Reprojection:', proj2);
var modisResolution2 = modisyr.projection().nominalScale();
print('MODIS resolution:', modisResolution2);


texts:
I have a python code that is generating coordinates one by one. Is there a way to display them directly in ArcMap as they are generated? I am aware of the insert cursor but the lock on the table prevents python from writing into it.

Ok, now that I used the Python window in ArcMap new error has occured. Could it be because the "czk" is a list?

Here is the code:

czk=[413960.0, 47130.0]
cursor = arcpy.da.InsertCursor("C:/GIS/AlJeba/Tocke/Tocke.shp", ("SHAPE@XY"))
for row in czk:
    cursor.insertRow(row)
del cursor

And here is the error:

Runtime error
Traceback (most recent call last):
  File "", line 333, in
TypeError: argument must be sequence of values


texts: Deriving touch points layer of line-polygon intersection in QGIS or PostGIS?

texts:
Assume you have polygon and a line that goes though this polygon.
Now I would like to create a seperate layer which stores the point geometries of where the line intersects the outline of the polygon, so where it enters or exists the polygon.
Is this possible in QGIS 3.x or PostGIS?


texts: I have just installed the Azimuth and Distance plugin for QGIS and I have been testing it for use with our COGO workload. However when using the Anti-Clockwise ArcDirection method for a curve to the left for example draws funky.  A curve to the right (clockwise works just fine). Has anyone had good luck COGO'ing parcels with legal descriptions including multiple curves left and right? You can see the overlapping hexagon is supposed to be the curve to the left. I am confused.


texts: What should i do if the Repository tab "Status" was indicating "unavailable".? (using proxy server)

texts:
In the plug in Repository tab. The "Status" was indicating "unavailable". I tried "hawkens" suggestion about the HTTPproxy thing but it didn't work (the network that I am using is a proxy server).

Any suggestion please? I'm trying hard to get the proxy host number and the port (is it the number on the host and port?)


texts:
We have this large map of United States on our website and we are trying to see it well on mobile. As the map is too large we lose some functionalities on mobile view like dragging.

Is there any possibility to have dragging capability for mobile using cartoDB?

I've changed this property to the code:

    map = new L.Map('cartodb-map', {
        center: [36, -96],
        zoom: 4,
        zoomControl: false,
    })

    // Enable drag handlers.
    map.dragging.enabled();

But is working just for desktop.


texts: What part of the Census Block is the County?

texts:
I have a list of census blocks for California and I am trying to get the county for each row.

If I have a census block of 60070004024003, which 2 digits make up the county FIPS?


texts:
I have this filter:

var filter_unset = new OpenLayers.Filter.Comparison({
        type: OpenLayers.Filter.Comparison.EQUAL_TO,
        property: attribute,
        value: 0
 });

Now, I want the filter to apply on null values as well, something like:

        value: 0 or NULL

Can I do this?


texts:
I have a number of shapefiles in different CRSs (mostly WGS84 lat/lon) that I'd like to transform into a common projection (likely Albers Equal Area Conic, but I may ask for help on choosing in another question once my problem gets better-defined).

I spent a few months doing spatial stats stuff in R, but it was 5 years ago.  For the life of me, I cannot remember how to transform an sp object (e.g. SpatialPolygonsDataFrame) from one projection to another.

Example code:

P4S.latlon <- CRS("+proj=longlat +datum=WGS84")
hrr.shp <- readShapePoly("HRR_Bdry"), verbose=TRUE, proj4string=P4S.latlon)
# Shapefile available at
#   http://www.dartmouthatlas.org/downloads/geography/hrr_bdry.zip
#   but you must rename all the filenames to have the same
#   capitalization for it to work in R

Now I have a SpatialPolygonsDataFrame with appropriate projection information, but I'd like to transform it to the desired projection.  I recall there being a somewhat unintuitively-named function for this, but I can't remember what it is.

Note that I do not want just to change the CRS but to change the coordinates to match ("reproject", "transform", etc.).

Edit

Excluding AK/HI which are annoyingly placed in Mexico for this shapefile:

library(taRifx.geo)
hrr.shp <-
  subset(hrr.shp, !(grepl( "AK-" , hrr.shp@data$HRRCITY ) |
                                     grepl( "HI-" , hrr.shp@data$HRRCITY )) )
proj4string(hrr.shp) <- P4S.latlon


texts:
GeoServer UI (http://localhost:8080/geoserver/web/) provides the means to configure the security of layergroups by granting read/write access to various roles. I would assume that the UI is simply a front-end on the REST API and that everything which can be achieved through the UI, can also be achieved via REST API directly.
I've tried a GET request on http://localhost:8080/geoserver/rest/security/acl/layergroups but received an error 404. This is to be expected, for according to the following API references, there is no endpoint for layergroups ACL:

https://docs.geoserver.org/latest/en/user/rest/api/accesscontrol.html
https://docs.geoserver.org/latest/en/api/#1.0.0/security.yaml

It would be nice if there were a separate endpoint for layergroups ACL. Apparently, after enabling read/write access of myrole to mylayergroup via GeoServer UI, the layergroups ACL settings show up in a GET request to the layers ACL:
GET http://localhost:8080/geoserver/rest/security/acl/layers
{'*.*.r': '*', '*.*.w': 'GROUP_ADMIN,ADMIN', 'mylayergroup.r': 'myrole', 'mylayergroup.w': 'myrole'}

However, GeoServer seems to expect layers ACL rules in the following format:
{'myworkspace.mylayer.r': 'myrole', 'myworkspace.mylayer.w': 'myrole'}

I tried a POST request with JSON payload in the same format as the UI seemed to have used, and got an error 422:
POST http://localhost:8080/geoserver/rest/security/acl/layers
JSON: {'mylayergroup.r': 'myrole', 'mylayergroup.w': 'myrole'}

    <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>
    <title>Error 422 Invalid rule mylayergroup.r, the expected format is workspace.layer.mode=role1,role2,...</title>

HTTP ERROR 422

Problem accessing /geoserver/rest/security/acl/layers. Reason:

    Invalid rule mylayergroup.r, the expected format is workspace.layer.mode=role1,role2,...

      Powered by Jetty:// 9.4.18.v20190429

How to configure the layergroups ACL via REST API?


texts: Difference between geometry(MultiPolygon,4326) and geometry(Geometry,4326)

texts:
I had imported a shapefile to PostGIS which has a column data type of geometry(MultiPolygon,4326).
I have used the following SQL statement before to find a point in the polygon of given column.
select precinctso
from public."precinct" b
where ST_Within(ST_GeomFromText(concat('POINT(',103.84263486000195,' ',1.2739257004756506,')'),4326),b.geom::GEOMETRY)
limit 1;

I have tried using ST_Contains:
select precinctso
from public."precinct" b
where ST_Contains(b.geom, ST_GeomFromText(concat('POINT(',103.84263486000195,' ',1.2739257004756506,')'), 4326))
limit 1;

But both statements return empty.
These statements have been used for a similar query where the data type of the column has been geometry(Geometry,4326).
I want to know whether the Geometry type MultiPolygon is any different when using ST_Within or ST_Contains.


texts:
I found a similar question, but the answer solve the problem indirectly. So, here I reshape the question.
Say, in the following situation, how am I going to use "ELSE IF" statement?
CASE
  WHEN  "land use" != 'residential' THEN 'other land uses'
  ELSEIF  "area" < 1000 THEN 'small lot'
  ELSE 'large lot'
END


texts: Is it possible to export layers scale-based visibility to a webmap?

texts: Is it possible to make or convert a QGIS map into an interactive webpage map? For example, I have layers that appear at certain zoom levels, is there a way to have that interaction on the map for HTML/Webpage?


texts: How to check for PostGIS Function Support Matrix for new custom EPSG imported from EPSG site

texts:
I need a support for CRS that was not available in PostGIS 2.5 so I looked up into the EPSG site they had simple easy straight forward insert query to add EPSG so I did that now how do I check or know that all the functions that PostGIS have will run with same accuracy as they use to work or may be I should avoid use of some specific function like given here , and all reprojection will work as required?

EPSG site reference https://spatialreference.org/ref/epsg/4979/

INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 94979, 'epsg', 4979, '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ', 'GEOGCS["WGS 84",DATUM["World Geodetic System 1984",SPHEROID["WGS 84",6378137.0,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0.0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.017453292519943295],AXIS["Geodetic latitude",NORTH],AXIS["Geodetic longitude",EAST],AXIS["Ellipsoidal height",UP],AUTHORITY["EPSG","4979"]]');


texts: How to move a layer position on the layers order

texts:
I am looking for a way of moving a layer to the top of the layer order once another layer is loaded. Using the following code extracted from this answer I get the error: AttributeError: 'QgsLayerTreeMapCanvasBridge' object has no attribute 'customLayerOrder'

    shapeindex = self.iface.mapCanvas().currentLayer()
    bridge = self.iface.QgslayerTreeCanvasBridge()
    order = bridge.customLayerOrder()
    order.inster(0, order.pop(order.index(shapeindex.id())))
    bridge.setCustomLayerOrder(order)

I looked into pyqgis3 documentation but it doesn't show more than the function name: QgsLayerTreeMapCanvasBridge


texts: Which parameter data type is used in the Integrate geoprocessing tool for ranking?

texts:
I'm trying to figure out which parameter data type is used in, for example, the Integrate geoprocessing tool.

It looks like this:

I'd like to use it in a geoprocessing tool I'm making, but can't figure out which data type it is.


texts: Connect / Disconnect GPS device via PyQGIS

texts:
In QGIS I would like to establish a connection to my GPS device programmatically (not via button in the GPS control panel)

I found out how to get available ports with QgsGPSDetector i.e.

QgsGPSDetector.availablePorts()
[(u'localhost:2947:', u'Lokaler GPSD'), (u'\\\\.\\COM4', 'COM4:')]

and how to get information via QgsGPSConnection to a connected device (cp. Accessing GPS from QGIS 2.14.1/Python 2.7 /Windows10).

Via

connectionRegistry = QgsGPSConnectionRegistry().instance()
connectionList = connectionRegistry.connectionList()

i can do

connectionList[0].close() rp. .connect()

But if no device is initially connected via the GPS control panel, the .connectionList() return nothing, so I'm not shure how to connect / disconnect to a GPS device.

ADD/EDIT

tried

con = QgsGpsdConnection('localhost',2947,'local gpsd')
# or con = QgsGpsdConnection('localhost',2947,''), rsp.
con.connect()
>>> True
con.status()
>>> 1
cr = QgsGPSConnectionRegistry().instance()
cr.connectionList()
>>> []

I'd like to try this 'COM4', but the port parameter is an integer.

Connecting via GPS Information panel does not work with gpsd, only Autodetect and serial (COM4:):

Do I have to construct a QgsGPSConnection that i can

.connect()

and

.close()

?


texts: How to Create 'Update Triggers' to Track Field Values Changes When Data is Versioned?

texts:
How can I create a trigger to capture when a field is changed from one value to another if my data is versioned? I am not able to run an UPDATE trigger and grab ':new' and ':old' values because with versioning the change in the field value is just an addition to the A#### table. I would like to be able to compare the old and new value of this field when it is changed to make sure they are not the same.

Using:

ArcGIS Pro
Versioned Database (I am working with and have access to a child checked out from DEFAULT as well as a grandchild (child of child of DEFAULT)
Oracle w/ SDE

I was able to create a trigger that captures the creation of a new feature in a specified feature class and stores it in a history_tracking_table (the start of tracking a record). It grabs a unique ID, the initial tracking state (which is NULL), the user who made the change, and the change date are captured as well. I have done this and it's functional but maybe not pretty. See code below:

create or replace
trigger insert_obs_history_2
AFTER INSERT
ON A261570
FOR EACH ROW

DECLARE
username varchar(20);
prev_Match_Recs number;

BEGIN
  SELECT COUNT(*) INTO prev_Match_Recs
  FROM HISTORY_TRACKING_STATE
  WHERE HISTORY_ACTION = 'INSERTED' AND
  HISTORY_REPORTID = :new.O_REPORTID;

IF prev_Match_Recs = 0 AND
:new.O_REPORTID IS NOT NULL AND
:new.SUBMISSION_STATE_ID IS NULL AND
:old.SUBMISSION_STATE_ID IS NULL
  THEN
    SELECT USER INTO username FROM dual;

    INSERT INTO HISTORY_TRACKING_STATE VALUES (
    :new.O_REPORTID, :new.SUBMISSION_STATE_ID,
    username, SYSTIMESTAMP, 'INSERTED' );
  END IF;
END;

How can I do something similar to capture updates to the SUBMISSION_STATE_ID field as it changes over the records lifetimes while working in a versioned environment where 'updates' are inserts into the delta tables?


texts: Address Listed on Open Data but cannot geo-locate

texts:
I am researching graffiti incidents in New York City.

One of the addresses entered is Frederick Avenue And West 179Th Street,Bronx,New York

Cannot geocode, i.e. find the latitude and longitude of this

I tried googling this, but even google maps cannot find.

How do I deduce where this can be.

For instance, there is an address entry, South East Corner Of Amsterdam Avenue And Cross Bronx Expressway,Manhattan,New York, and I can approximate it to the exit of Cross Bronx Expressway that takes one to Amsterdam Avenue, i.e. Exit 1.

How to resolve this issue?


texts: How to import a non-spatial Table into QGIS?

texts:
Can I import tabular data in Qgis, without any information about X and Y only the table?
I need to join them after importing them.


texts: Extract the polygon layer vertices and delete the unwanted ones in QGIS

texts:
I need to extract the vertices coordinates of the four corners of each rectangle of a layer made up of 729 polygons (= map sheets). I began using  Vector/Geometry Tools/Extract Vertices function to create a vertices layer but each rectangle contains 7 vertices - the first and the last are superimposed, I suppose - so I want to eliminate the median vertices and one of the two (the first or the last) vertices superimposed. Moreover the 4 vertices of each rectangle coincide with the vertices of the adjacent rectangles having the same coordinates but different attributes related to the different map sheets.

Now I do not know how to delete the unwanted vertices (I can not use the Processing / Vector General / Delete Duplicate Geometries function because it would also delete the vertices of adjacent rectangles).
How to solve this problem?


texts: Is there a way to programatically trigger popups using the operationallayer config?

texts: I have a widget that queries a group of features, and I'd like to be able to trigger a popup that looks and behaves the same as the one that gets generated when the user clicks on the map. Is there a call to do this in the flexviewer itself?


texts: Reduce the resolution of satellite imagery data for use in MATLAB

texts:
When I am using the 4-band multi-spectral data for NDVI calculation, after successful completion of runtime, it shows some warning message like "Loading 17%..."

How should I rectify it? Should I do any filtering for reducing the resolution effect before processing?


texts: How to download large subsets of USDA NAIP imagery from Box.com?

texts:
The US Department of Agriculture’s Geospatial Data Gateway states, ‘As of April 21, 2017 the NAIP datasets are only available through the "NAIP Download" option on the home page and are no longer be available through the Gateway ordering process.’ NAIP (the National Agricultural Imagery Program) is a dataset of satellite imagery of the continental US, and the link they provide lands me on this Box.com shared link.

Box.com limits downloads to 15 GB, so the only way I can download all of 2017 North Dakota data is to manually download files one-by-one, or select a few of them that are hopefully less than 15 GB and download as several batches.

I’ve asked on Box’s help forum and thought to also ask the GIS gurus here: is there any way to automatically download this data, either through Box or some other way?

Non-solutions I've tried:

WebDAV doesn't follow shared links that I save in my Box.com (I signed up as a free user).
FTP is not available for free users.
The Box REST API doesn't recognize folder_ids that I get from the URL (i.e., at the list of 2017 data, the folder for nd or North Dakota is 43484391906, but 43484391906 is not a valid folder_id according to the API).
Disabling JavaScript gives a list of files and directories but no download link works (looks like Box.com is only usable by JavaScript-enabled clients).
This makes me wonder if using the Box CLI will also turn out to be a non-solution in getting this data automatically.

Is it possible that the data is not accessible through the REST API because it's hosted on an on-premise appliance? (The nrcs.app prefix in the URLs makes me think this.)


texts:
I am using QGIS 2.18 and plugin qgis2web to produce a web map.
In QGIS I have the layer layout I want, with a particular image sitting above all others. This layer-order is changed when I attempt to generate the web map (see image).
What am I doing wrong?
Is there a way to dictate layer order in the web map?


texts: I would like to create numerous line segments from a set of GeoJot points. Is there a way I should format the data so I can run a script; and is there a script that can accomplish this? I have coordinates (geotagged photos) for Intersection "A" and "B" and would like to link them with a line that would then have the data form from the point (deplicate data form for "A" and "B" except one says "begin" and the other "end"; I don't care about this field, so could use either point). I have done an entire community this way and it is tedious to create new segments by manually inserting a line between the points, moving the data form to the line and then deleting the points.


texts: is there any analogue of the ArcGIS "topo to raster" tool in QGIS? I need to creates surfaces from shp for some cases from points and for some from lines and this instrument is very convenient.


texts: append X,Y,Z and intensity values into a las file

texts:
I want to append XYZI values row after row into a las file.How to do this?
I'm using python laspy library to write XYZI points at once, but i want to append row after row using a loop.

The points are in SQLite db I need to generate a .LAS file(v1.3 and v1.4).I get memory error as the db grows large.

i'm doing this

import laspy

header = laspy.header.Header()
outfile = laspy.file.File("output.las", mode="w", header=header)
outfile.X = read column1 from db
outfile.Y = read coulmn2 from db
outfile.Z = read column3 from db

outfile.close()


texts:
I am doing some basic work with G Sheets and Sheetsee.js library.

I have it working as far as I get my coordinates showing up on a map. I'm using the demo-map file as my starter, I just changed the URL key and added my own column names to optionsJSON function.

My problem is, none of my column data is being called into the pop up template. I cannot find out what I'm doing wrong.

Here's my code for the script part:


texts: ST_Buffer in meters from 4326 point table

texts:
I would like to create 100 meters circle buffers in a UTM projection from a point table whose SRID is EPSG:4326 and keeping the id column.
I wonder what is the straightest way. Is that possible to get it in one query after having created a buffer table with the URM SRID? Like:
INSERT INTO buffers SELECT id_points, ST_Buffer (geom, 100 , 'quad_segs=8') AS geom from points

but this would require I guess to add ST_transform somewhere and ST_Buffer does not like to involve two different SRID.


texts: Missing band in mapping function using Google Earth Engine

texts:
I built a function to add a a new band (b1_1) to each image of a image collection and I mapped it. It does seem to work, at least when printing the image collection and inspecting it in the console.

But when I trying to map on the "new collection" another function for charting  the new band (b1_1) against some other pre-existing band, the new one doesn't exist anymore.

//mask and add delta DEM band
var add_band = function(image){
var glacier_img = image.select(['b1']);
var glacier_clip = glacier_img.clip(shp);
var dem_glacier = DEM_delta.clip(shp);

var vel_range = glacier_clip.lte(ee.Number(700)).and(glacier_clip.gte(ee.Number(-700)));
var dem_glacier_mask = dem_glacier.updateMask(vel_range);

var final_img = image.updateMask(vel_range).addBands(dem_glacier_mask);
return final_img;

};

//add a DEM band  to each image of the images collection
var imgCollection = imgCollection.map(add_band);
print(imgCollection);

//chart
var graph = function(image){
var id = image.id;
var image = ee.Image(image.id);
//print(id);

var result = image.reduceRegion(ee.Reducer.toList(),shp, 200);
var y = ee.Array(result.get('b1'));
var x = ee.Array(result.get('b1_1'));
var chart = ui.Chart.array.values(y, 0, x)
      .setSeriesNames(['b1'])
      .setOptions({
        title: 'Velocity vs. Delta_DEM ',
        hAxis: {'title': 'Delta_DEM'},
        vAxis: {'title': 'Velocity'},
        pointSize: 3,
        });

  // Print the chart.
  return print(chart);
  };

 imgCollection.evaluate(function(imgCollection) {  // use map on client-side
 imgCollection.features.map(graph);
 });


texts: How to add application/json as GetFeatureInfo format in a WMS with QGIS Server?

texts:
I have a qgis-server instance.
In my GetCapabilities request I can see these formats available for the GetFeatureInfo:

<GetFeatureInfo>
   <Format>text/plain</Format>
   <Format>text/html</Format>
   <Format>text/xml</Format>
   <Format>application/vnd.ogc.gml</Format>
   <Format>application/vnd.ogc.gml/3.1.1</Format>
[...]

I would like to add the application/json format.
What is the approach?

Can I use the qgis plugins?


texts: FME: How to remove duplicate points in linestrings using a tolerance on geometry

texts:
I have a set of linestrings which are contaminated by points being almost duplicate.

Example:

  X------XX--------X-X-----X----XX-------X

I would like to remove the duplicates using a tolerance on the geometry, leaving me with a set of clean linestrings:

  X------X---------X-X-----X-----X-------X

I know its probably a simple question, but being new to FME i could use a simple educative hint or answer :)


texts: QGIS DBmanager: error creating the SQL layer

texts:
i've converted an OSM file from Geofabrik to an SQLite database and loaded it into the QGIS DBManager (QGIS 3.12). I've executed the following SQL query and got a nice table as result.
SELECT *
FROM points
WHERE place = "city"

Yet when I try to load the layer into the canvas using the Load button (and with the settings seen in the screenshot), I receive the following vague error:

There was an error creating the SQL layer, please check the logs for further information.

There's nothing to be seen in the logs.
Does anyone know how to fix this?


texts:
I am trying to update a popup for a polygon that I have drawn with Leaflet.Draw. Now I understand that to set a popup to a layer I just created, it's all fine. I can also get the id of that layer I just created using the L.stamp(my_layer).

Let's say I create 3 layers: A, B and C. Each of them has their stamp id stored in my database. Now let's say at some point I want to bind  a popup to layer C. How do I access that layer C when I click on it and bind the popup to C only?


texts: QGIS attribute table: round values but don't change the data

texts:
In my attribute table of a vector layer I have the following field with the decimal data.

I would like that QGIS shows me only the first 2 decimal positions in the Attribute table, but I do NOT want QGIS to actually change the data.
Just to show me the data a bit nicer.
Is that possible? Is there a function I'm not seeing?
I'm using QGIS 3.4


texts: Free alternative to "Depression evaluation" tool from ArcHydroTools

texts:
Is there any free program or package that has the same function as Depression evaluation tool in ArcHydro tool from ArcGIS?

I am currently working on extracting depressions from DEMs as well as their respective drainage area. Depression analysis function does exactly what I want, which is defining the depressions in DEM and the drainage area of each depression. However, I no longer have access to ArcGIS.

I know that I can extract the depressions using the stochastic depression analysis tool from Whitebox tools. Also, I am familiar with QGIS, but I was unable to find the same thing in either QGIS or Whitebox.


texts: Getting the majority value for a number of rasters in a conditional statement

texts:
I am trying to execute a conditional statement using Python and ArcGIS 10.1 that looks something like this:

outRas = Con( MAJORITY (raster1, raster2, raster3, raster4) > 110, 1, 0)

The idea is that I want to check if the most common value for the current pixel among the provided rasters is greater than 110, and if so, assign a value of 1 to the output raster. Of course, there is no such function called MAJORITY, however I'm wondering if there is something equivalent, or if anyone knows of any other way to achieve this.


texts: Exporting toolbox from ArcGIS 10.5 (Python 2.7) to ArcGIS 9.3 (Python 2.5)?

texts:
I built an interface using ArcGIS 10.5.0 with Python ver. 2.7, I want to export this toolbox to an other PCs that using ArcGIS 9.3 with Python ver. 2.5.

I did it using Toolbox save as but the toolbox seems to be empty.


texts: "invalid 'times' value" error in stack function

texts:
I would like to stack multiple GeoTIFF file. My code looks like this:
var1 <- list.files(path="D:/#This is the access path", pattern=".tif$", full.names=TRUE)

var1

var2 <- stack(var1)

and here is the error:

Error in rep.int(factor(names(x), unique(names(x))), lengths(x)) :
invalid 'times' value


texts: problem with graduated styling in qgis - not rendering one feature

texts:
I am trying to style a polygon layer as graduated according to a specific field of the polygon shapefile in QGIS.

I selected equal interval or equal count and then sort but when I clicked OK one feature is not rendered. Its value equals the upper limit of the upper class. This occurred with the other types of styling like pretty jerks too. What is the problem?


texts: Map Algebra - Grass GIS (r.mapcalc) merge two raster with NoData value

texts:
I have a slope raster and a road raster. The black pixels in the road raster are NoData values and the pink pixels have a value of 1. They don't cover exactly the same area.

I want to create a new raster based on the slope raster, except for pixels where there is a road pixel. If there is a road pixel, I want to assign a value of zero to the slope raster.

I am not quite sure how the formula needs to look like.
I tried the following, put it didn't give me the expected results:

r.mapcalc "output = ((Slope+(Road*null())))"

Any ideas?


texts:
i'm trying to export image from GEE to my drive and I have gotten this error messege:
Error: Exported bands must have compatible data types; found inconsistent types: Byte and Float64.

I have read here that it may have to do with the bands, but the thing is that I don't really need all the bands, I want to download it in RGB .

this is the code I have :

// Load a Landsat 8 -atmospheric corrected
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
//surface reflectance
    .filterBounds(table)
    .filterDate('2015-03-01', '2018-03-31');

// Get the number of images.
var count = collection.size();
print('Count: ', count);

// make a composite (mean image) of the images
var meanImage = collection.mean().clip(table);
print(meanImage,'mean image');

//make image for mesopotamia only
//var mesopotamia=meanImage.clip(table);

//visualization
var visParams = {
  bands: ['B4', 'B3', 'B2'],
  min: 0,
  max: 3000,
  gamma: 1.4,
};

Map.addLayer(meanImage,visParams,'Landset8');
//Map.addLayer(mesopotamia,visParams,'Mesopotamia');

var toexport = meanImage.visualize(visParams).addBands(meanImage)

  // Export the image, specifying scale and region.
Export.image.toDrive({
  image: toexport,
  description: 'MesopotamiaLandset8',
  scale:30,
  maxPixels:201406799,
  region: table
});