texts: Fetching single value without Cursor in ArcObjects?

texts:
I used the below function to get the single value using ArcObjects.

private string getApproverAlias(string checklistId, bool isDeleted)
{
    string approverAlias = "";
    var masterTable = featureWorkspace.OpenTable("TableOne");
    IQueryFilter masterTableFilter = new QueryFilterClass();
    masterTableFilter.SubFields = "ApproverAlias";
    masterTableFilter.WhereClause = "ChecklistId = '" + checklistId + "' AND IsDeleted = " + isDeleted;
    ICursor cursor = masterTable.Search(masterTableFilter, true);
    IRow row = null;
    while ((row = cursor.NextRow()) != null)
    {
        approverAlias = row.get_Value(row.Fields.FindField("ApproverAlias")).ToString();
    }
    return approverAlias;
}

I know from the above method the output is single value, but I'm using Cursor and Row to fetch the result.

Is it any other way to get the result without using the Cursor and Row in ArcObjects?

Something similar to ExecuteScalar in ADO.net like the below:

SqlConnection connection = new SqlConnection("Data Source=ServerName;Initial Catalog=DatabaseName;User ID=UserName;Password=Password");
connection.Open();
SqlCommand command = new SqlCommand("SELECT ApproverAlias FROM TableOne WHERE ChecklistId = '" + checklistId + "' AND IsDeleted = " + isDeleted, connection
string approverAlias = command.ExecuteScalar();
connection.Close();


texts:
I am trying to create a loop so that each layer is produced in its own map.

This is my code:

import arcpy
path = "C:\Users\Et\Documents\ArcGIS\Default.gdb"
mxd=arcpy.mapping.MapDocument(path +"\WinterGritting.mxd")
df=arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
lyrList =arcpy.mapping.ListLayers(mxd,"",df)[0].name
elements=arcpy.mapping.ListLayoutElements (mxd)
for eachlyr in lyrList:
    arcpy.SelectLayerByAttribute_management(eachlyr,'NEW_SELECTION')
    df.zoomToSelectedFeatures()
    arcpy.SelectLayerByAttribute_management(eachlyr,'CLEAR_SELECTION')
    mxd.title=eachlyr.name
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToJPEG(mxd, path + "\Routes" + eachlyr.name + "jpg")

del mxd

The Error I get is:

  Traceback (most recent call last):   File
  "C:/Users/Et/Documents/Python", line 9, in
      arcpy.SelectLayerByAttribute_management(eachlyr,'NEW_SELECTION')   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 6494, in
  SelectLayerByAttribute
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000732: Layer Name or Table View: Dataset L does not exist or is
  not supported Failed to execute (SelectLayerByAttribute).


texts:
What is currently the highest value of OSM_ID in openstreetmapping?

What is the highest possible value of OSM_ID? 2^31 (2.147.483.648)?

And here I add some more text or my question is rejected because it is too short...


texts: Updating records in a table using GetCount and an Update Cursor

texts:
What I have is a buffer which is then used to select a number of point features which intersect the polygon. I would like to take that number and add the value to a table. The row where the value will need to be located corresponds with the name of the original feature class the buffer was created from.

import arcpy
from arcpy import env
import os
from os.path import join

arcpy.env.workspace = r'MyWorkspace'
updateTable = r'table'
shortField = "FieldName"
arcpy.MakeFeatureLayer_management("FieldName", "fieldname_layer")
fcs = arcpy.arcpy.ListFeatureClasses("*_Value", "POLYGON")
for fc in fcs:
    select = arcpy.SelectLayerByLocation_management("fieldname_layer", "INTERSECT", short)
    result = int(arcpy.GetCount_management(select).getOutput(0))
    print short.replace('Buffer_','') + " " + str(result)

So far the best I can come up with is printing the file name less the Buffer_ portion which precedes the name I need. I cannot wrap my head around how to go from there to updating only relevant records using the update cursor.


texts: Restricting WMS GetFeatureInfo properties in OpenLayers

texts:
I am working on a WMS-based OpenLayers project, which is pulling in data from the BGS WMS Servers. We're pulling the 1:50k bedrock layer from this WMS:

https://map.bgs.ac.uk/arcgis/services/BGS_Detailed_Geology/MapServer/WMSServer?

In OpenLayers 3 I am trying to query the WMS and display the attributes in an iFrame beneath the OpenLayers map. Using the following code, I am able to get the full metadata to display beneath the map.

map.on('singleclick', function(evt) {
            document.getElementById('info').innerHTML = '';
            var viewResolution = /** @type {number} */ (view.getResolution());
            var url = wmsSource.getGetFeatureInfoUrl(
                evt.coordinate, viewResolution, 'EPSG:3857',
                {'INFO_FORMAT': 'text/html'});
            if (url) {
              document.getElementById('info').innerHTML =
                  '<iframe seamless src="' + url + '"></iframe>';
            }
          });

The thing is, this pulls in ALL of the metadata and I would like to restrict it to just one or two fields, such as LEX, LEX_RCS_D or MAX_TIME_D (see field names here)

I'm not a user of JavaScript, so I'm having difficulties modifying the code.

After some research online, I came across the propertyName parameter, which I've tried implementing, although I'm not even sure this is what I'm looking for.

{'INFO_FORMAT': 'text/html'}, {propertyName: 'LEX'}

Given the way the BGS WMS is configured, is it even possible to limit the attributes which are pulled out on the GetFeatureInfo request?


texts: Reoccurring multipart / slivers / topology errors using file geodatabase?

texts:
Recently I had to hand in a FGDB to state authorities. I did all topology checks and corrections as defined in their specifications (see list below). My collegue did a cross check on a the same FGDB. We are a hundred percent sure that all errors and topology issues were fixed. We handed in the data and it came back with a list of errors.

We thought this cannot be true and ran another check on our data (exactly the same checks on the data we handed in). We couldnt believe our eyes, the errors mentioned by the authorities were reproducable.

How can this be possible? What have we possibly done wrong? Could it be a server problem? Or are the somewhat strange tolerances given in the original FGDB responsible for such problems? Has anyone experienced similar issues?

Required topology checks:

no multipart features --> we used dissolve / explode multiparts (editing session) / multipart to singlepart tools
no features smaller than 20 m2 --> we used the eliminate features tool
FC1 must be inside FC2, FC2 must be inside FC3, etc. --> topology checks

Our work environment:

ESRI ArcGIS 10.3, Advanced License
FGDB on a server
cluster tolerance required for topology errors: 0.01 m
however, the empty FGDB sent to us by the authorities that needed to be filled with data had the following tolerances set:


texts: Can we read/write MXD Document using GeoTools

texts: Is it possible we can read/write mxd document only using geotools.


texts:
I'm trying to create a table with spatial column in oracle and insert data into it. I'm using ArcSDE and Oracle. The version of Oracle is 11.2.0, the version of ArcSDE is 10.0.

My create statement is:

CREATE TABLE sensitive_areas (area_id integer, name varchar(128),
area_size float, type varchar(10), zone sde.st_geometry);

And insert statement is:

INSERT INTO SENSITIVE_AREAS (area_id, name, area_size, type, zone)
VALUES (1, 'Summerhill Elementary School', 67920.64, 'school', sde.ST_PolyFromText('polygon ((52 28,58 28,58 23,52 23,52 28))', 0));

I've already added st_shaplib.dll and libst_raster_ora.dll to the environment variable EXTPROC_DLLS.

My tnsname.ora:

    LISTENER_ORCL =
  (ADDRESS = (PROTOCOL = TCP)(HOST = localhost)(PORT = 1521))

ORACLR_CONNECTION_DATA =
  (DESCRIPTION =
    (ADDRESS_LIST =
      (ADDRESS = (PROTOCOL = IPC)(KEY = EXTPROC1521))
    )
    (CONNECT_DATA =
      (SID = CLRExtProc)
      (PRESENTATION = RO)
    )
  )

ORCL =
  (DESCRIPTION =
    (ADDRESS = (PROTOCOL = TCP)(HOST = localhost)(PORT = 1521))
    (CONNECT_DATA =
      (SERVER = DEDICATED)
      (SERVICE_NAME = orcl)
    )
  )

My listener.ora:

SID_LIST_LISTENER =
  (SID_LIST =
    (SID_DESC =
      (SID_NAME = CLRExtProc)
      (ORACLE_HOME = D:\app\rainman\product\11.2.0\dbhome_3)
      (PROGRAM = extproc)
      (ENVS = "EXTPROC_DLLS=ONLY:D:\app\rainman\product\11.2.0\dbhome_3\bin\oraclr11.dll;D:\Program Files\ArcGIS\ArcSDE\ora11gexe\bin\st_shapelib.dll;D:\Program Files\ArcGIS\ArcSDE\ora11gexe\bin\libst_raster_ora.dll")
    )
  )

LISTENER =
  (DESCRIPTION_LIST =
    (DESCRIPTION =
      (ADDRESS = (PROTOCOL = IPC)(KEY = EXTPROC1521))
      (ADDRESS = (PROTOCOL = TCP)(HOST = localhost)(PORT = 1521))
    )
  )

ADR_BASE_LISTENER = D:\app\rainman

However, I got this error: ORA-28595: Extproc agent : Invalid DLL Path. What should I do?

Thanks

USER_LIBRARIES entry for ST_SHAPELIB is also correct:
ST_SHAPELIB D:\Program Files\ArcGIS\ArcSDE\ora11gexe\bin\st_shapelib.dll    Y   VALID


texts: Getting path of project, or layer file in PyQGIS?

texts:
I am generating new raster files from vector layers on a project. I would like to save them in the same directory as the project or layer files.

How do I find this path using PyQGIS?


texts: Is there a way to break up an image file into separate polygons based on colour?

texts: I have a set of land-use maps that I created years ago in ArcMap (see image attached), I need to do more work using them but I have since lost the shapefiles. Is there any way in ArcGIS that I can break this image into polygons based on the different colours of the image? Or do I need to redraw them all by hand?


texts: PostGIS calculate distance to every point in table

texts:
I would like to calculate distance from point made by putting coordinate manually to every point in table. I made that code:

SELECT a.name
FROM ne_110m_populated_places as a
CROSS JOIN (SELECT ST_Point(20, 50)::geography AS ref_geog) As r
WHERE ST_DWithin(a.geom, ref_geog, 10000000)
ORDER BY ST_Distance(a.geom, ref_geog);

And it is working fine but I don't know how to add distance to output table as a column with value.


texts: Different image quality for different zoom levels (OpenLayers/GeoServer)?

texts: Original data that raster layers are made of are of a high quality. At a small scale (1:100) this is useful. However at a large scale (1:10 000) it results in ~300KB PNG images which are actually not needed to be of that good of a quality. I am looking for an approach to have cached raster layers where one group of zoom levels (lets say range from 11 to 20) images would have different (lower quality) than it should based on original high quality data.


texts:
I want to draw a dashed line on straight line

i tried but (Example)

* {
   stroke: yellow, red;
  :nth-stroke(1) {
     width:4,
     dasharry 5 2
   };
   :nth-stroke(2) {
     width:5;
   }
 }


texts: How to put zipcodes data onto a web map?

texts:
This is to create a cellphone coverage map. I'm being given a list of zipcodes and need to create a detailed static cellphone map (plus, ideally, an interactive zoomable map) of cellphone coverage areas.

I'm totally new to GIS systems, or any mapping API other than Google Maps. Can you suggest how I would get started?


texts: How should I handle CRS properly after the major change in PROJ library?

texts:
As explained here, the PROJ library has changed a lot and one should stop using proj.4 strings. To me it is not clear, what I should do instead. This answer explains, how I could deal with the problem without any changes in one concrete use case, but to me it seems, this is more a hack than a good idea on a long term. (At the moment, this is the only thing, I can do. The only alternative is not to update important packages.)
Is there any advice for the following process:
For example, I download data (e.g. boundaries of cantons) as geojson-files from map.geo.admin.ch.
As we don't want to do the process each time, we import the data and save it as RData-file:
boundaries <- "my_boundaries.geojson"
data <- geojsonio::geojson_read(boundaries, what = "sp")
# Warnmeldung:
# In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
# Discarded datum Unknown based on Bessel 1841 ellipsoid in CRS definition,
# but +towgs84= values preserved
data2 <- sp::spTransform(
  data,
  sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
save(data2, "boundaries.RData")

When I analyse data according to polygons, I do the following:
boundaries <- load("boundaries.RData")
data_coords <- data.frame( # from some source mydata
  lat = mydata$lat,
  lng = mydata$lng, stringsAsFactors = FALSE)

spdf <- sp::SpatialPointsDataFrame(
  coords = data_coords,
  data = data_coords,
  proj4string = raster::crs(boundaries))
df_in <- sp::over(spdf , boundaries , fn = NULL)

So my question is:

What is the proper process for GIS-data after the major change in PROJ library?
What does this mean in our concrete case above?
I also found a nice example with plots, but it also uses proj4string. What should I do with this example?

(Unfortunately, I am not a GIS-expert)

Using CRS throws an error, if used as suggested in the answer by Roger Bivand below:
data2 <- sp::spTransform(
  data,
  sp::CRS(SRS_string = "OGC:CRS84"))
# Fehler in h(simpleError(msg, call)) :
#   Fehler bei der Auswertung des Argumentes 'CRSobj' bei der Methodenauswahl
#   für Funktion 'spTransform': Fehlender Wert, wo TRUE/FALSE nötig ist

This is working well:
spdf <- sp::SpatialPointsDataFrame(
  coords = data_coords,
  data = data_coords)
slot(spdf, "proj4string") <- slot(boundaries, "proj4string")


texts: Accessing smb share on Trimble YUMA with python

texts:
I need to write a script that would import a .shp on the YUMA to a gdb on the script/host pc for manipulation with arcpy. I can't seem to find a smb python module. So Im not entirely sure on how to do this.

Have I just not found it? surly its came up before.


texts: Extracting raster values with polylines in QGIS

texts:
I have a polyline shapefile that represent a road an a raster grid with slope values.
Is it possible to interpolate the shape file and the grid to have the average slope of the road?
I know that I can interpolate grid with points and polygons. But with lines?


texts: How to set a new raster band to a QgsPalettedRasterRenderer object?

texts:
Is there a way to set the current raster band to an existing QgsPalettedRasterRenderer object? Other renderer have specific methods for that (e.g. setRedBand, setGreenBand, setBlueBand, setGrayBand, etc.).

For QgsPalettedRasterRenderer I haven't found anything like that. I guess it is done in a different manner (?).

renderer = layer.renderer()
assert isinstance(renderer, QgsPalettedRasterRenderer)

# Getting the current raster band is working:
print(renderer.band())

# But how to set a new raster band? A setBand method does not exist.
renderer.setBand(band) # not working


texts:
How can I apply italic fonts for labels and print composer in QGIS 3.4?

The option seems to be disabled.


texts: Zonal statistics as table in Model Builder changes zone field names

texts:
I am using Model Builder to calculate the mean NDVI for some polygons and then I want to combine those tables to the attribute table of those polygons. For whatever reason, the attributes in the Zone Field in the zonal statistics as table tool keep changing so the model fails because either 1) whenever I try to join the table based on the attribute the model stops because the join field is not recognized because the zone field changed or 2) the zone field changed and "is not recognized".
I select FID_vegeta (see image 3), hit ok, then when I run the model and it fails and when I open the tool back up and the zone field now says name of polygon layer.FID_vegeta which I never selected because it was not an option.


texts: Detect cattle on high resolution satellite images

texts:
I have no deep experience in remote sensing. I would like to know if it might be feasible to detect cattle (i.e. a single cow or bull) from a high resolution satellite image e.g. Worldview 3  with a spatial resolution of 0.31m. I was thinking about visual interpretation. Big animals have a diameter of about 1.80 meters and their colour (white brown) is clearly distinguishable from the ground below (green).

The aim is just to be able to identify if there is one or more animals (generally the herd is bigger than just a single animal).

This question is rather speculative since I don't think anybody has probably done anything similar before but maybe you could indicate if you think I should give it a try with a sample image...


texts: pyqgis: change transparency for cells with value 0

texts:
I have a Raster file in which some cells have the value 0 and I would like to use pyqgis make them 100% transparent. I guess I have to use the QgsRasterTransparency class but I have no idea how.

I have a QgsRasterLayer object called rlayer and I can make a QgsRasterTransparency obejct, but I don't know how to apply the transparency on the rlayer object and how to define that the transparency only applies to those cells with value 0.


texts:
I use Google maps in my application but because this service has limitation (2500 request per day for Geocoding) I want to use OSM instead of Google maps.

But OSM doesn't cover Iran as well as Google maps. My question is: Can I use Google maps tiles in OSM for Iran?

You can compare Google maps tiles with OSM in Iran on the following pictures: First is OSM and second is Google maps.


texts: Splitting line into segments of different length with ArcGIS Desktop?

texts:
I have a river network and I have to split each river into transects of different length. The distances where the lines have to be split are in an excel sheet.

For example I have to split
the river 1 at the distance of 200, 400, 1000 m;
the river 2 at the distance of 100, 500, 1300 m;
and so on...

I would like to do it with ArcGIS Desktop (I have standard license not advanced).

After some attempts I could locate the points where I want to cut my line. To do that I treated my lines as routes, I saved my excel sheet as dBASE to use it as event table. Then with the function 'make route events layer' I created a shapefile with the point where I have to cut the lines.
Now the problem is that I should probably use the tool 'split line at point', but I dont have the advanced license.

Is it possible to do it with another tool using a lower level license?


texts: Create and immediately update field using Field Calculator

texts: You can create fields and once created, populate them with data using the Field Calculator. But can you do both at the same time using an expression? Or must you do it in 2 steps?


texts: Doing difference between two layers in PyQGIS?

texts:
I'm developing a plugin where I would like to get some features from one layer (a), do a difference on them with the entirety of a second layer (b), and put the results into a layer I've created in memory (c). If at all possible, I'd also like to copy a specific field from (a) to (c) when doing the difference.

To do this, I've got the following code set up:

layer_A = None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "A":
        layer_A = lyr
        break

it = layer_A.getFeatures( QgsFeatureRequest().setFilterExpression ( u'"CAT" = 27') )
layer_A.setSelectedFeatures( [ f.id() for f in it ] )

box = layer_A.boundingBoxOfSelected()
iface.mapCanvas().setExtent(box)
iface.mapCanvas().refresh()

layer_B = None
for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
    if lyr.name() == "B":
        layer_B = lyr
        break

l1 = QgsVectorLayer("Polygon?field=CAT:integer", "C", "memory")
QgsMapLayerRegistry.instance().addMapLayer(l1)

As you can see, we select various features in A based on a query.

I'd like to perform a difference between these selected features, and the entirety of B, which is unfiltered and to put the resulting shapes into C (with the CAT, if at all possible, but not necessary).


texts: Is there a way to quickly paste multiple columns of data into an attribute table in QGIS?

texts: So I have a shapefile with 133 rows of data, and I want to add 4 extra fields to it. My problem is that after I have created the new fields in the attribute table, I have no way to quickly paste the data for these columns from my spreadsheet. Can someone show me how to do this?


texts: Can I use external TMS service as a store in Geoserver?

texts:
I have GeoServer hosted on my server and I want to connect a store so that I can serve layers from there. I am aware that I can connect external WMS source to geoserver, but not sure how to connect TMS.
e.g. of TMS Service
https://sit-geonetwork.geoplatform.info/geonetwork/srv/eng/catalog.search#/metadata/ca63648a-2330-5df7-b644-e2a53832acda


texts: Using python elif syntax in Field Calculator to reclassify a range of values from a different field and populate a new field

texts:
I'm trying to use the Field Calculator with a python elif statement to reclassify values from a different field and populate a new field, and while I'm not getting an error message of any kind, I'm also not getting any sort of result.  In other words, the tool runs and it says that it succeeded but nothing has changed in the column that is being calculated.

In codeblock example below the column that I want to use the elif statement on is "pa_gNyr_ha" and the new column that I want to populate is "N_load_bin".

def Reclass(Nload):
  if (Nload >= 0 and Nload <= 10000):
    return 5000
  elif (Nload > 10000 and Nload <= 24000):
    return 17000
  elif (Nload > 24000 and Nload <= 40000):
    return 37000
  elif (Nload > 40000 and Nload <= 58000):
    return 50000
  elif (Nload > 58000 and Nload <= 74000):
    return 66000
  elif (Nload > 74000 and Nload <= 90000):
    return 82000
  elif (Nload > 90000 and Nload <= 104000):
    return 97000
  elif (Nload > 104000 and Nload <= 118000):
    return 111000
  elif (Nload > 118000 and Nload <= 130000):
    return 124000
  elif (Nload > 130000 and Nload <= 142000):
    return 136000
  elif (Nload > 142000 and Nload <= 157000):
    return 149500
  elif (Nload > 157000 and Nload <= 173000):
    return 165000
  elif (Nload > 173000 and Nload <= 189000):
    return 181000
  elif (Nload > 189000 and Nload <= 210000):
    return 194500
  elif (Nload > 210000 and Nload <= 233000):
    return 221500
  elif (Nload > 233000 and Nload <= 261000):
    return 242000
  elif (Nload > 261000 and Nload <= 299000):
    return 280000
  elif (Nload > 299000 and Nload <= 356000):
    return 327500
  elif (Nload > 356000 and Nload <= 456000):
    return 406000
  elif (Nload > 456000):
    return 501500

And then the expression is

Reclass(!pa_gNyr_ha!)

Essentially what I'm trying to do is create a conditional statement that will look at a value in the "pa_gNyr_ha" field and return a value to the "N_load_bin" field based on the range it falls into.

Am I at least close?


texts: Rename a empty layer that has been added to a group layer

texts:
Originally I was trying to add a empty layer that was stored on disk. That contained no data. I created the empty layer , by making a group layer within my mxd and then saving it onto a location on my computer

The process was detailed here

Create empty group layer within group layer

My original question was "once my layer has been added to my group layer . How do I rename it "

At the time I was struggling to understand how to access empty layer I had just added. I was getting bogged down in trying to understand the list feature of the group and its sub groups. Thinking I need to access the Empty layer first by going into its top level group. But In this case I did not need to list the layers because I was all ready able to access the layer from the code that had loaded into my Group

So all that is needed to add to the code , as supplied on the previous post is

mxd = arcpy.mapping.MapDocument(r"C:\Temp\test.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
groupLayer = arcpy.mapping.Layer(r"C:\Temp\Group.lyr")

Inside here you can then rename the layer you another name

groupLayer.name = 'NEW_NAME'

arcpy.mapping.AddLayer(df, groupLayer, "BOTTOM")
targetGroupLayer = arcpy.mapping.ListLayers(mxd, "Group", df)[0]
addLayer = arcpy.mapping.Layer(r"C:\Temp\Group2.lyr")
arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, addLayer, "BOTTOM")
mxd.save()
del mxd,df


texts: Why is R writeRaster function with GTiff format so slow?

texts:
I've downloaded a countries shapefile from Natural Earth, and used it to clip GLC2000 raster using Brazil's border as the mask (in R).

setwd('~/R/GLC2000')
library(rgdal)
library(raster)
glc <- raster('glc2000_v1_1')
pa <- readOGR('../shp/ne_10m_admin_0_countries','ne_10m_admin_0_countries')
p <- grep('Brazil',pa$ADMIN)
p1 <- crop(glc,extent(pa[p,]))
glcP <- mask(p1,pa[p,]) # dividing the mask in crop+mask makes it much faster
dim(glcP)
plot(glcP)
writeRaster(glcP,as.character(pa@data[p,'ISO_A3']),'GTiff')

dim(glcP) gives 4369 5056    1, which is big, but not absurdly big.

plot(glcP) gives the image below, which proves it worked as expected.

So why doesn't writeRaster complete its task, even hours after being issued?

I want to to this for every country in the world, that's why I'm trying to do it in R. It's easier for me than a script in, say, QGIS.


texts: QGIS Map Tip can't display photos on the network computer

texts:
I'm using Map Tip to display image but if photo store in my computer it display ok ,example:

img src="C:\test.jpg"

But if photo store on a computer in my local network it can't display ,example:

img src = " \\computer1\Folder1\photo.jpg "

I don't know how to solve it.


texts: Cluster processing resources needed for 50 Gb DEM?

texts:
For "renting purposes", I have to specify the computing resources I need from the computing center at my University. I have very little idea of what to ask for.

Considering the information below may you suggest something?

This is the computational cluster:
"A ten blade server with each blade consisting of 2 Quad Core Intel Xeon 2.33 Ghz CPUs with 16 GB of memory for a total of 80 CPUs and 160 GB of memory. Accounts on a 2 TB disk system. Running 64-bit Linux"

With lossless compression, the original data is one 50gb tiff file.
At some point in time i shall be working with several (up to ten) files 50gb in size.

I would use GDAL, Python scripting and maybe C++ scripting. If allowed, I would use softwares like Grass GIS and Saga GIS. Also, R language with spatial libraries.
I will be deriving the usual terrain parameters, trying to apply object-oriented algorithms for the extraction of specific features (landforms) and using statistical analysis for descriptive and modeling purposes.


texts: Key error when trying to add a new attribute value to feature - PyQGIS3

texts:
I get this error KeyError: 'FieldOne' when I try to add a new attribute value to my feature.
This is the code that I use to create a new Layer:

v_layer = QgsVectorLayer('LineString?crs=epsg:4326', 'layerName', 'memory')
pr = v_layer.dataProvider()
pr.addAttributes([QgsField('FieldOne',  QVariant.String)])
pr.addAttributes([QgsField('FieldTwo',  QVariant.String)])
v_layer.updateFields()
QgsProject.instance().addMapLayers([v_layer])

After creating the new layer, I need to fill it with new features. Each feature has his attributes. Using this code that give me the error KeyError: 'FieldOne'.

def createNewFeature(fieldone,fieldtwo):
  seg = QgsFeature()
  seg.setGeometry(QgsGeometry.fromPolyline([start_point, end_point]))
  pr.addFeatures( [ seg ] )
  v_layer.startEditing()
  seg['FieldOne']=fieldone
  seg['FieldTwo']=fieldtwo
  v_layer.updateFeature(seg)
  v_layer.commitChanges()


texts: How do I load multiple shapefiles into the same existing feature class?

texts:
I have a hundreds of shapefiles and, based on the filename, I need to load it into one of three feature classes (points, polylines, or polygons). My code works until it tries to load the second file of a particular type. Then, I get an error "FC_NAME already exists". Yes, it already exists, that's exactly what I want. Now load the shapefile into the feature class, please!

How do I load a shapefile into an existing feature class?

Here's the code, trimmed down to exclude logic that gets to this point:

result = arcpy.FeatureClassToFeatureClass_conversion(shp, outputLocation, outputPointsFC)


texts: Map layer disappears when I change projection

texts: I am using qgis 2.16.1. When I change the projection (using Layer>Set CRS of layer(s), the map disappears. I have attempted to recover the map by changing back to the original projection with no success. How can I change the projection and maintain a visible map?


texts: Digitize objects in a LiDAR point cloud?

texts:
I am working with LiDAR point cloud data of road networks. The goal is to digitize all visible objects like hydrants, streets, sidewalks.
My input data is available in LAS or ASCII format.
My first thought was to use ArcScene and the 3D Analyst to use the points as backdrop data to digitize the features. But I don't have any license yet. See: Displaying LAS Dataset in ArcScene.
I am wondering if there are other software which you can use to digitize point clouds. Many of them can do automatic classification, but that is not what I want. The determining factor is not "free software" here.

I've installed the free Trial of ArcGIS Pro and tested the 3D Editor on my LAS dataset. I was able to draw features within the point cloud. Nonetheless, the usability and performance was pretty bad for a data set with 50,000,000 points. Seems like ArcGIS is not the best tool to draw 3D objects with geographic information.


texts: Improving accuracy of address returned by Google Maps using reverse geocoding?

texts:
We are using mobile application to track the location of the users. we have done testing of the application and found that the accuracy of the address is not correct most of the time. while testing we have also found that multiple users on the same location getting different addresses from google which is sometimes 1km away from their actual locations.

we are using reverse geocoding to get the address of the user. mobile application sends the latitude and longitude to server. server makes the reverse geocoding api call to the google and fetch the address from the google.

I want to know what we can do to improve the accuracy of the results returned by google using reverse geocoding?


texts: ArcMap sometimes ignores labeling expression

texts:
I'm using this python script as an expression for labeling in ArcMap:
def FindLabel ( [volkorrekt] ):
 if [volkorrekt] >= '5':
  return "<FNT size = '20'>"+  str( [volkorrekt] )+"</FNT>"
 elif [volkorrekt] >= '2':
   return  "<FNT size = '14'>"+  str( [volkorrekt] )+"</FNT>"
 elif [volkorrekt] >= '1':
  return "<FNT size = '8'>"+  str( [volkorrekt] )+"</FNT>"
 else:
  return ""

I want ArcMap to only label the feature if the value volkorrekt is bigger than 1, this is working. I also want ArcMap to change the font size according to my if-statement (font size 20 if value >= 5, 14 if >= 2 and 8 if >= 1) but this only works for some but not for all features.
This is my current result:

Examples for wrong labels: The "19,487..." is displayed very small though it's one of the biggest values, the "9,552..." however is displayed much to large.
How can I get ArcMap to label the features in the way I described?


texts:
I have a script with an bail-out function like so:

def die():
    from sys import exit
    exit()

Occasionally in the script I check a condition and exit if necessary. This works great in 10.1, but unfortunately our server is still at 10. I had to modify the function as follows to work (per this post):

def die():
    from sys import exit
    try:
        exit()
    except SystemExit:
        pass

The problem is that it gets to the except, passes, and keeps on trucking right to the end of the script. Is there any way to actually make the script stop?


texts: Learning to think spatially rather than just follow point and click tutorials

texts:
The college I work for is considering revamping their intro to GIS course. Traditionally, they've used the ESRI suite but are now open to using other options. There are two factors responsible for the desire to make a change:

The text book provided by ESRI. The instructors feel that the point and click tutorials don't teach the students how to think about GIS but rather how to use ESRI software. They'd like to have something that's more than a manual and that won't be out of date by the time they arrive.
Traditionally graduates could expect jobs from either the Territorial Government or Parks Canada (both of which use ArcGIS), but with recent gov't cut backs these two are no longer the main destination for graduates. After consulting with industry representatives the college has been told that they should be training student to "think spatially" rather how to use just one type of software.

Are there any resources the SE.GIS community would consider helpful for developing GIS thought and instruction (rather than tutorials)? And if you were going to commit to teaching with one program/suite which do you think would be most useful to learn such that the skills learned would be most broadly applicable to GIS software in general?


texts: How do you store/organize/manage your own data?

texts:
Lately I was wondering how one could or should store his/her own data.
As I studied GIS and I am working on different private projects using GIS data, a lot of data have accumulated on my local hard disk. Up to now, I have been collecting all data in a simple folder structure. Now I was thinking if that is the correct way or if it would be better to store the data in a geodatabase (Postgres with Postgis) or on a server (Geoserver) on my localhost.
The main problem is, that I have created folders that store the data thematically (forest, water, basemaps, ...), but I would like to give my data multiple tags, so that I can filter them not only thematically, but also based on the study area, year, license and so on.
The data are vector and raster data.
How do you store, organize and manage your data?


texts: Openlayers Reprojection from 4326 to 3395 is wrong?

texts:
I'm using proj4js definitions downloaded directly from epsg.io webpage.

When trying to reproject a small extent I get this result:

ol.proj.transformExtent([-180, -90, -179, -89], 'EPSG:4326', 'EPSG:3395')

Result: [-20037508.342789244, -30198185.16987704, -3.12413936106985, -1.5707963267948966]

Isn't this absolutely wrong?


texts: Edit listbox and table value from a dialog box

texts:
Good day! I have a code that import my table on a listbox and able to add on it.

i was stock on some codes when i double click a selected value on my listbox to open another dialogbox to edit (the selected value) it as well as on my table.

Heres my code...

Include "Mapbasic.def"

Declare Sub Main
Declare Sub AddItemDialog
Declare Sub AlterMainDialog
declare sub editlistDialog
declare sub Nodelist(byval sId as string)
dim listboxValues() as string
declare sub listbox_edit

'---------------------------------------------------------------------------
Sub Main

   Dialog Title "Main Dialog" Width 150 Height 80 calling editlistDialog
      Control ListBox       Position 13, 7      Width 100 Height 50 ID 1 Title from variable listboxValues calling listbox_edit
      Control Button Title "Add Item" ID 2 Calling AddItemDialog
      Control OkButton
      Control CancelButton

End Sub
'---------------------------------------------------------------------------
Sub AddItemDialog

   Dim strAddItem as String

   Dialog Title "Add Item Dialog"
      Control EditText Width 100 Into strAddItem
      Control OkButton Title "Add"
      Control CancelButton

   If CommandInfo(CMD_INFO_DLG_OK) Then
      Dim i as Integer
      i = Ubound(listboxValues)
      Redim listboxValues( i + 1)
      i = Ubound(listboxValues)
      listboxValues(i) = strAddItem

      insert into subn(lname)
      values(strAddItem)

      Call AlterMainDialog
   End If

End Sub
'---------------------------------------------------------------------------
Sub AlterMainDialog

   Alter Control 1 Title from variable listboxValues

End Sub
'---------------------------------------------------------------------------

sub editlistDialog
    call Nodelist("5350")
    alter control 1 title from variable listboxValues

End Sub
'---------------------------------------------------------------------------
sub Nodelist(byval sId as string)
    dim i as integer
    dim nRows as integer

    select * from subn where site=sId into lstValues noselect
    nRows=tableinfo("lstValues",tab_info_nrows)
    redim listboxValues(nrows)
    for i = 1 to tableinfo("lstValues",tab_info_nrows)
        fetch rec i from lstValues
        listboxValues(i)=lstValues.lname
    Next
End Sub
'---------------------------------------------------------------------------
Sub listbox_edit
    Dim i as string
    If CommandInfo(CMD_INFO_DLG_DBL) = TRUE Then
    ' ...then the user DOUBLE-clicked.
    ' see which list item the user clicked on.
        i = ReadControlValue(1)
        print i
        'Should try to look the table and edit list

    End If

End Sub

br,

Ivan


texts:
I am following a tutorial on color-relief hillshades (found here: https://ieqgis.wordpress.com/2015/04/04/create-great-looking-topographic-maps-in-qgis-2/ )

The post is dated 2015-04-04, so - not that old.

The author says "superior colour ramps can be downloaded using Etienne Tourigny’s Color Ramp Manager (Plugins – Manage and Install Plugins)."

I cannot find this in the Plugin manager. It must be a current plugin, because there is a page for it at https://plugins.qgis.org/plugins/colorrampmanager/ .

How can I find & install this plugin, please?

I am using QGIS 2.12 Lyon on a Mac running OS X 10.10.5.


texts: Deducing start and end of season information using r.seasons in GRASS GIS at 'Field' level

texts:
I would like to ask one question on using r.seasons in GRASS GIS. The information page on this function is available at https://grass.osgeo.org/grass78/manuals/addons/r.seasons.html
The problem that I am facing is with regards to deducing seasonal information of Sugarcane Crop for a particular area between 2016 and 2020. I am using Sentinel-2 imagery and have already got the decadal maps (10-day maps). I am doing the analysis at the 'Field' level with demarcated boundaries, whereas the output of the r.seasons function is at the pixel level.
Using the zonal statistics function combined with majority mode, I am able to deduce the number of seasons in every field but when I am trying to fetch the start and end of the seasons, there is an overlap because of the different seasonal properties of different pixels. This overlap got ignored when the majority mode was applied but to carry out further analysis, I require the following information at field level to forecast the crop yield:

The number of seasons for every field.
The start and end decadal of every season in every field.

How can I solve this confusion?


texts: Exporting RGB values from raster to file that can be processed in Excel using ArcGIS Desktop?

texts: I have a 3 band raster file with RGB values (the Black Marble night time imagery from NASA). I want to get these values out and into excel so that I can plot them and hopefully see some sort of clustering for some values. I do however not know how to get the RGB values exported to a file format that works in Excel.


texts:
I have map images of my local ward (Areas geo-referenced and drawn using QGIS and Adobe illustrator. Now these maps are in .png format (Raster images). It's not a tile linked to some mapserver such as OSM or Google Maps, but plain vector images which have routes "drawn" out.

I want to be able to link these routes to independent pages (Each route has a set of data stored in a mysql db). I want people to be able to hover over these routes on the vector image with a highlight effect for the route (routes can be curved or straight). And when they click on a particular route, it opens a popup which gives details about the route pulled from the db.

I understand that the PHP GD library can be used for this purpose, but how do i highlight/link routes using this library? If anyone can point me to specific functions from this library, that would help.

Help would be greatly appreciated.


texts:
I am trying to map the density of buildings in flood-prone areas. In order to do that, I created a grid of equal-sized squares. I want to intersect the layer of buildings with the grid so I can determine the density of at risk buildings in each square.
Problem is, the intersection doesn't seem to work.

On ArcGis : I get an feature class but only some of the buildings have been intesected, others are still astride two squares.

On QGis : I get an empty output.

Does anyone have an idea of what the problem could be ? Am I using the wrong tool or is there a problem with the data I am using ? ArGis also refuses to extract the data from selection layers I created from the buildings... I guess it could be damaged, but I don't know how to fix it.


texts: Export a QGSVectorLayer as KML format with style

texts:
I have a vector layer with geometry type PointGeometry and styled with QgsGraduatedSymbolRendererV2 renderer.

I have tried:

QgsVectorFileWriter.writeAsVectorFormat(layer,"kml_example.kml","utf-8",None,"KML")

But it only exports points i.e. without colors.

How can I export this vector layer as KML output with colors?


texts:
I want to keep updated my fresh install of Nominatim 3.1.0. I have installed North America maps from this (http://download.geofabrik.de/north-america-latest.osm.pbf) file.

According to the article (http://nominatim.org/release-docs/latest/admin/Import-and-Update/) I can use 'CONST_Replication_Url' in settings/local.php for daily country diffs update but I can't find one update file like in the article's example ('https://download.geofabrik.de/europe/ireland-and-northern-ireland-updates'). For North America I go to (http://download.geofabrik.de/north-america/) and see us-midwest-updates, us-northeast-updates, us-pacific-updates, us-south-updates and us-west-updates folders...

Anybody knows where I can get one update file?

Because now I use this command ./utils/update.php --import-osmosis-all and the update process works all time...like a live update and I don't know actually this update works for only installed maps or for the whole world.


texts:
I have a problem with QGIS. When using the contour tool I get the following error:

  The process failed to start. Either the invoked program is missing, or
  you may have insufficient permissions to invoke the program.


texts: "The Fields collection did not contain an expected geometry field" (ArcGIS 10)

texts:
I am writing a fairly simple application to build simple polylines.

I got most of the application working for the most part, but this part is just baffling me to no end and there is nothing I can find about this problem, so maybe I am alone on this.

This is the code:

        IFeatureClassDescription fcDescrip = new FeatureClassDescription() as IFeatureClassDescription;
        IObjectClassDescription objDescrip = (IObjectClassDescription)fcDescrip;

        IFieldsEdit fields = (IFieldsEdit)objDescrip.RequiredFields;

        IFieldEdit geoField = (IFieldEdit)fields.get_Field(fields.FindField(fcDescrip.ShapeFieldName));

        IGeometryDefEdit geoDefField = (IGeometryDefEdit)geoField.GeometryDef;

        geoDefField.GeometryType_2 = esriGeometryType.esriGeometryPolyline;
        geoDefField.SpatialReference_2 = spatialReference;
        geoDefField.HasM_2 = true;
        geoDefField.GridCount_2 = 1;
        geoDefField.set_GridSize(0, 0);

        for (int i = 0; i < fields.FieldCount; i++)
        {
            Console.WriteLine(fields.get_Field(i).Name);
            if (fields.get_Field(i).Name == "SHAPE")
            {
                Console.WriteLine(fields.get_Field(i).GeometryDef.GeometryType);
            }
        }

        UID uidCL = new UID();
        uidCL.Generate();
        UID uidEXT = new UID();
        uidEXT.Generate();

        IFeatureClass featureClass = featureWorkspace.CreateFeatureClass(m_name, (IFields)fields, uidCL, uidEXT, esriFeatureType.esriFTSimple, string.Empty, string.Empty);
        IFeature feature = featureClass.CreateFeature();

It tells me that it doesn't contain a geometry field... even though it clearly does because it tells me what shape it has in there.

I tried using add_Field, but it does the same thing (oddly) if I do. It clearly updates 'fields' since I had it output the type before and after and it changes from polygon to polyline.

I looked at this thread and I tried defining everything, but it still gives the error.

Any ideas?


texts:
I am trying to simplify several polygons part of one shapefile using Python and OGR. However, I am stuck with how to use the simplify function as mentioned in the documentation:

http://www.gdal.org/classOGRGeometry.html#a0f4d7948332d9efd6548e1cf87bb6c8f

I tried this, for example:

data = ogr.Simplify(inData)

But then I get an AttributeError:

  'module' object has no attribute 'Simplify'

Any tips?


texts:
I am trying to put road No in the shield with tilemill .mss. Here is the code I used.

#road {
    shield-file:  url(C:\wamp\www\shield.png);
    shield-name: [ROAD_NO];
    shield-face-name:"DejaVu Sans Bold";
    shield-size:10;
    shield-fill: white;
    shield-spacing:50;
    shield-line-spacing:5;
    shield-type: png;
        }

Its return two errors in the shield-file and shield-name fields.

What wrong I have done here?


texts:
I'm testing the ability to use SQLLite with ArcGIS. I ran the python script Esri said, "sqlite_database_path = 'C:/sql_lite/CI.sqlite' and  arcpy.gp.CreateSQLiteDatabase(sqlite_database_path, "ST_GEOMETRY") it looks like it worked, ArcCatalog sees the database.

What I found is you have the options of using spatialLite, ST_Geometry, or the geopackage. I was thinking the ST_Geometry would let you use sqllite directly.

I saw comments where people copy/pasted data into sqllite but in ArcCatalog I didn't have the option to paste after I copied. I tried the ArcPy Copy_management, and import, both screamed "the application is not licensed to create or modify schema for this type of data"

I also found Esri said:
You cannot create a geodatabase in SQLite.
Feature services are not supported with SQLite.

So how are people putting GIS Featureclasses into sqllite from ArcMap?


texts: Can I change the GDAL/OGR library that QGIS uses?

texts: I have been adding a new format to OGR, all seems well and OGRInfo works. Now I would like to test my code out using QGIS. Is there an easy way to switch which GDAL lib QGIS uses on a windows box? I've read Tim's blog post but would rather not go down the recompiling from source route if possible.


texts: Adding new vertex using QGIS Vertex editor

texts:
I'm having trouble adding new vertex using the QGIS3 Vertex Editor - I'm using
QGIS version 3.4.1-Madeira.  I tried following the instructions provided in answer to Using QGIS vertex editor?.

The instructions stated:

  QGIS 3.0 - Two options to add a new vertex:

  Double-click on a segment to add a new vertex. The new vertex
  starts out selected and stuck to the cursor. Click a third time to
  place the vertex.
  When hovering over a segment, there is an X in the middle of the
  segment. Click on the x to add a new vertex. The new vertex starts out
  selected and stuck to the cursor. Click a third time to place the
  vertex.

I've tried both of these with the following results:

double clicking on a segment seems to select and move an existing vertex,

and

hovering over a segment does not result in an x in the middle of the segment that I can click and move.

I have tried both options 1 and 2 above with snapping enabled and disabled and this doesn't seem to make any difference.

Any suggestions as to what I'm doing wrong?


texts: ArcObjects: ExportOperationClass causing error and not exporting

texts:
I'm trying to export each row in a shapefile, but I get an error I haven't seen before (see below). I feel like it might be the labyrinth of interfaces that I have to cast between that may be causing the error. I must be getting close though because I get an export feature class progress bar (see below) which then gets overridden by the error.

Code:

try
        {
            IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
            IMap map = mxdoc.FocusMap;

            for (int i = 0; i < map.LayerCount; i++)
            {
                ILayer layer = map.Layer[i];
                if (layer.Name == "lga_polygon")
                {
                    if (layer is IFeatureLayer)
                    {
                        IFeatureLayer flayer = layer as IFeatureLayer;
                        IFeatureClass fc = flayer.FeatureClass;

                        IDataset pDataset = flayer as IDataset;
                        IDatasetName pDatasetName = pDataset.FullName as IDatasetName;

                        string lgaName = "LGA_NAME";

                        IFields fields = fc.Fields;
                        int lgaNameFieldIndex = fields.FindField(lgaName);

                        IField lgaNameField = fields.Field[lgaNameFieldIndex];

                        IFeatureCursor lgaCursor = fc.Search(null, true);

                        IFeature lga = lgaCursor.NextFeature();

                        while (lga != null)
                        {
                            string nameOfLGA = lga.Value[lgaNameFieldIndex].ToString();

                            IDatasetName dsName = new FeatureDatasetNameClass();
                            dsName.Name = nameOfLGA;
                            IFeatureClassName fcName = dsName as IFeatureClassName;

                            IWorkspaceName pWorkspaceName = new WorkspaceNameClass();
                            pWorkspaceName.PathName = "C:\\Users\\username\\Documents\\ArcGIS\\Default.gdb";
                            pWorkspaceName.WorkspaceFactoryProgID = "esriDataSourcesFile.FileGDBWorkspaceFactory";

                            IQueryFilter qFilt = new QueryFilterClass();
                            qFilt.WhereClause = "LGA_NAME = '" + nameOfLGA + "'";

                            IGeographicCoordinateSystem pGCS = new GeographicCoordinateSystemClass();
                            SpatialReferenceEnvironment spatRefEnv = new SpatialReferenceEnvironmentClass();
                            ISpatialReference spatRef;
                            pGCS =  spatRefEnv.CreateGeographicCoordinateSystem(Convert.ToInt16(esriSRGeoCSType.esriSRGeoCS_GDA1994));

                            spatRef = pGCS as ISpatialReference;

                            IGeometryDef geoDef = new GeometryDefClass();
                            IGeometryDefEdit geoDefEdit = geoDef as IGeometryDefEdit;
                            geoDefEdit.SpatialReference_2 = spatRef;
                            geoDefEdit.GeometryType_2 = fc.ShapeType;
                            IGeometryDef geoDef2 = geoDefEdit as IGeometryDef;

                            IExportOperation exportOp = new ExportOperationClass();
                            exportOp.ExportFeatureClass(pDatasetName, qFilt, null, geoDef2, fcName, ArcMap.Application.hWnd);

                            lga = lgaCursor.NextFeature();
                        }
                    }
                }
            }
        }
        catch (Exception ex)
        {
            MessageBox.Show(ex.Message + "\n\n" + ex.Source + "\n\n" + ex.StackTrace);
        }

Error:


texts: Local geolocation to FIPS county using Latitude and Longitude

texts: I have millions of lat,long pairs which I want to convert to FIPS county codes that match with the US census data. In the past I used the census block api but that is too slow for the amount of numbers that I have. Is there a better library / tool that I can load locally ? If not would be possible to use the TIGER data from the census website to create a service like that ? What knowledge is required (projections, minimum distances, ...) ?


texts: How to create random points on a line in QGIS

texts: Can I generate random points on a shapefile created from a digitized road? If so, how can I do that?


texts: Identifying line intersects within single polyline shapefile using QGIS?

texts:
I have a shapefile representing a road network, there is an issue with the network in that some roads have duplicated features in the shapefile (seperate digitisations with seperate attributes but for the same road).

My objective is to identify these duplicates.

As the attributes are different, I can't look in there for the duplicates so I was wondering if this could be done spatially?

I can't find a tool which allows me to check intersects within one file though, only between two seperate shapefiles.


texts: PostGIS query that selects feature if it goes outside boundary

texts: I am looking of a postGIS function that will select the highlighted line in the picture below. Using st_intersection I can clip the lines to the polygon, but i'm looking to select the one line that breaks out of the polygon. The trouble is the line originates inside the polygon so st intersects, st crosses, st intersection both return all lines not just the one that breaks out of the polygon.


texts: Checking if point is within multipolygon using ST_Contains?

texts:
Why is this query always returning false?

SELECT ST_CONTAINS(the_geom,(select ST_GeomFromText('POINT(26.8715 80.7799)'))) FROM blockgeom

Can I not check if a point is within a multipolygon using ST_Contains?


texts: i.landsat.toar not working and returning error

texts:
First time using grass. I need to convert DN to reflectance. Imagery -> Sat. images tools-> i.landsat.toar. But when I add raster name. I get error.
Im using landsat 8 images

WARNING: WARNING: ESUN evaluated from REFLECTANCE_MAXIMUM_BAND
Calculating...
ERROR: Unable to open header file for raster map <LC81890272015184LGN00_B11@abcdefgh1@

I have try this but it wasnt working.


texts: Change size to the squares of the confusion matrix

texts:
I have a problem with size in the 'plot_confusion_matrix', the squares of the confusion matrix appear cut off. I would like to solve this problem.

svc = SVC(kernel='linear',C=1,probability=True)
svc.fit(X_train, y_train)

model_evaluate(svc, (X_train, y_train), (X_test, y_test))


texts:
I am using NTS to generate a shapefile based on users' selection of data in a web interface. Each item selected is represented by a linestring and all relevant dbf data is exported nicely and the user ends up with a single shapefile that contains all of the lines selected.

I would like to be able to color code each line based on the data associated with each line but do not see anything in the NTS library that would allow me to do that. Is it possible to do this and if so, how? The solution has to support pretty much any GIS tool.


texts: Shapefile vs Arrays of coordinates to build Google Maps overlay?

texts:
I am very new to GIS, and my task is to overlay about 100 polygons on a map. I'm using Google Maps Javascript API v3.

I think I can define all my polygons as sets of lat/long coordinates in my database and programmatically loop through those sets of coordinates to create polygon overlays to put on the map.

Or, it looks like I can take an existing ESRI/ArcGIS shapefile and convert it to KML and use that for a custom overlay.  This seems more difficult to me and like it may be problematic if I need to tweak my polygon data if I don't like what's in the shapefile.  But I have no expertise with this, so I'm not sure.

At any rate, I don't know the advantages or disadvantages of either approach.  Or if there's another approach I should consider.

Any guidance or words of wisdom are greatly appreciated.


texts: Sextante - SAGA- Shape to grid question?

texts: I have a vektor shape (Index 1-5). After using the saga tool (shapes to grid) I'm getting a raster (values 1-5). The problem is that the NO DATA area is defined as 129. Is it possible to set 129 as no data? I cannot find the tool.


texts: how do you convert graphics to features and then associate those features with an existing layer?

texts:
Say for example, that you have water ties. These ties connect from a home to the main water line along the street. You want to draw these ties in ArcMap and give them attribute values regarding how long the ties are. However there is already an existing layer and many of the ties have already been drawn and connected to this layer. You are told to finish drawing the remaining lines and input the remaining measurements for those lines into the already established layer's attribute table.

After drawing a line from the Polygon (home) to the waterline, how do you relate it to the existing layer so that you can enter its length information?


texts: List of all fields and their domains for each table using SQL (table name, field name, domain name)

texts:
I'm trying to list all fields (in all tables in an Oracle geodatabase) and the coded value domains that they use.

+------------+------------+-------------+
| table_name | field_name | domain_name |
+------------+------------+-------------+
| table_1    | field_1    | domain_A    |
| table_1    | field_2    | null        |
| table_1    | field_3    | domain_B    |
| table_1    | field_4    | domain_B    |
+------------+------------+-------------+
| table_2    | field_1    | null        |
| table_2    | field_2    | domain_C    |
| table_2    | field_3    | domain_A    |
+------------+------------+-------------+

I've tried using the third example on this page: Query geodatabase system tables with SQL. It gets me the table names and the associated domains. But it's missing the field name:

SELECT reltypes.name AS type,
   origin_items.name AS "Origin Name", dest_items.name AS "Dest Name"
FROM sde.gdb_items origin_items,
   sde.gdb_itemrelationships relationships,
   sde.gdb_items dest_items,
   sde.gdb_itemrelationshiptypes reltypes
WHERE
   origin_items.UUID = relationships.originid AND
   dest_items.UUID = relationships.destid AND
   relationships.type = reltypes.UUID AND
   reltypes.name = 'DomainInDataset' AND
   origin_items.name LIKE 'PUBWORKS.%'
ORDER BY  origin_items.name

How can I make a list of all fields and their domains, for each table using SQL? Keep in mind that a domain can be used by more than one field, per table (or even by multiple tables).

Environment is Oracle 12c 10.3.1 Geodatabase.


texts: ArcGIS Modelbuilder: Merge features from same layer into one feature

texts:
How can I merge two features from the same layer into one feature when using the Model Builder?

Outside of the model builder, this can easily be done with the Edit Tool. How can I make the same process work in Modelbuilder?

Also, is there a way to automatically tell Modelbuilder to chose the Feature with the lower Value of Field-X to be the "Parent-feature"?

I work with ArcGIS 10.1
I also asked the question on Geonet: https://geonet.esri.com/message/608891


texts: How can I ask for location and display it on map on my website?

texts:
I have added a map to my website. Now I want to add a feature to it wherein when a user clicks on the add event button on the header a popup should appear which ask for event name and location and as soon as the user submits the information the location gets marked on the map with a style icon.

(I have added the code to show the popup with the relevant input form but I am not able to link the location with the map with a style icon containing the information).


texts: QField Sync v3.2.0 Photo Naming not available

texts:
I have QFieldSync 3.2.0 installed, but do not see the 'Photo Naming Tab' in Project Configuration that was provided in 3.1.8.
Has this been removed?


texts: How to pre-populate JavaScript forms in ArcGIS Web Application?

texts: I have an ArcGIS web application and would like to pre-populate order numbers.  Say if order 1 has been entered for a point on a map, when the next point is added, order 2 would populate automatically.  Of course these numbers are a series of letters and numbers, say abcd1234, next order would automatically populate abcd12345.  I am using SQL server to submit new features with feature access capabilities.  The web application is using the JavaScript API.  I have the code hosted on an IIS web server in order to access the actual JavaScript.


texts: Empty layer after import into PostGIS Database. Pyqgis

texts:
I'm facing troubles to import a layer into a PostGIS database with pyqgis in a script. Here is my code:

Nombre_BBDD='test'
host='localhost'
puerto='5432'
usuario='postgres'
password='xxxxxxx'
Nombre_Tabla='table'
Llave_primaria='id'
Capa_a_procesar='/xxxxxxx/shapefile.SHP'
Tipo_geometria='POLYGON'

capaImportar = QgsVectorLayer(Capa_a_procesar, 'capaImportar', 'ogr')
uri = """dbname='"""+Nombre_BBDD+"""' host='"""+host+"""'    port="""+puerto+""" user='"""+usuario+"""' password='"""+password+"""' key="""+Llave_primaria+""" type="""+Tipo_geometria+""" table=\"public\".\""""+Nombre_Tabla+"""\" (geom) sql="""
crs = QgsCoordinateReferenceSystem(25830)
error = QgsVectorLayerImport.importLayer(capaImportar, uri, "postgres", crs, False, True)

I was trying to do it how is recommended in this link.

It seems everything is right, but when I try to load the layer this is empty and the python console throws this error:

  Exception RuntimeError: 'wrapped C/C++ object of type PGVectorTable has been deleted' in bound method PGTableDataModel.del of db_manager.db_plugins.postgis.data_model.PGTableDataModel object at 0x7f3f24190d60 ignored

I'm working with qgis 2.18.10 in a Kubuntu 16.04


texts: Visually differentiating raster pixels in area of same value QGIS

texts:
Given a discrete band, and say an area that for example all has the value 5, is there a way in which I can highlight or differentiate where pixels start and stop?
This doesn't have to be some permanent differentiator among all pixels like some corner marking (although that might be helpful), even something that can just visually highlight the pixel I'm clicking on. I'm probably overlooking something simple but have been searching for a while and cannot find a solution.


texts: Labeling specific features in ArcGIS for Desktop?

texts:
In ArcGIS 10.3 for Desktop, I work with this point layer and I want to remove all the Sg & Sgs labels from it:

So I can get this result only:

All the labels are in one field named "field_name" (it is a string type field).
I can seperate the data into two fields or export the features with the numbers only.

I prefer to remove it with SQL Expression. I tried this one:

def FindLabel ([field_name]):
  if "Sg" IN [field_name] :
    pass
  return [field_name]

I am getting this error:

Error 0 no line 0.
SyntaxError: invalid syntax (<string>, line 2)


texts: Matching vectors in separate layers in QGIS. Overlay Question?

texts:
I am not exactly sure how to phrase this question since I do not know the term I need to use for this. I am planning driving routes to cover an entire state's road system. I already have a driving routes vector layer I have been creating and a state roads vector layer that includes all roads I need to cover. I need to know which roads have not been driven...Which state roads have not been traced over with my driving routes?

Overlay? Intersect? Join?


texts: gdal_contour is missing lines below a certain elevation

texts:
I have a very strange problem with gdal_contour. I have a very large GeoTif of Europe and would like to generate contour-lines for that region. Everything worked fine when I had lines every 25m but now I wanted to do it again and have lines every 10m.

It runs through without any error, but I am missing lines below 500-600m. Some lines until 500m are present, but many aren't. And below 500m I don't have any lines. I already ran it twice to just to make sure if there was a problem at first run, but second run had the same issues. Is there anything I am doing wrong here? Here's my command that I used.

gdal_contour -a height -f PostgreSQL -i 10.0 input "PG:host=localhost user=postgres password=postgres dbname=contours"

I need to have the output in a Postgres DB cause in a later process I would need to simplify some of the lines. So choosing a shapefile wouln't be an option.


texts: Unable to encode List data type as attribute when export FeatureCollection

texts:
var FC = ee.FeatureCollection([
  ee.Feature(ee.Geometry.Point(120, 0), {'instrumentMode': ['IW', 'EW']})
  ]);
print('FC', FC);
Export.table.toAsset(FC)

When I tried to export this table, I confronted the error message, saying that "Unable to encode value 'instrumentMode' of feature 0: invalid type". I am confuse why making a list as a atribute does not work?


texts: Applying ArcPy SpatialJoin_analysis in GeoPandas

texts:
I'm trying to implement the same solution as ArcPy's SpatialJoin_analysis in GeoPandas.
I have one shapefile which contains Polygons and Multipolygons. I need to perform below spatial operation of ArcPy using open-source GeoPandas library.
arcpy.SpatialJoin_analysis(
    target_features=spatial_join_in,
    join_features=spatial_join_in,
    out_feature_class=spatial_join_out,
    join_operation="JOIN_ONE_TO_ONE",
    join_type="KEEP_ALL",
    field_mapping="#",
    match_option="ARE_IDENTICAL_TO"
)

The Same shapefile is passed as both spatial_join_in files.
This operation also adds Join Count in each row, which is later used to filter the density areas.
How can we achieve similar solution using GeoPandas?


texts: Using ArcMap to build a Terrain DEM that looks exactly the same as a TIN built from the same point shapefile

texts: I am attempting to replicate DEMs previously constructed as TINs in the Terrain format. Building a Terrain and a TIN out of the same point shapefile produces two maps that look extremely similar, but when I run the results through the Surface Volume the Terrain and TIN give slightly different volume results. The volume of the Terrain also changes depending on whether Z Tolerance or Window Size is used as a pyramiding method even though the Surface Volume tool is analyzing the full resolution Terrain data. Is there a specific set of parameters that will produce a Terrain that gives identical volume results to a TIN or are the two methods inherently different in some way?


texts: OpenLayers Vector polygon and jQuery Droppable

texts:
How can I use jQuery Droppable to drag and drop an Icon over an openlayers polygon vector and select the polygon while over on hover? I have managed create a create a popup when over a div element, but I can't figure out how to do the same for a polygon inside the "map" div. Hope that makes sense.

Here is a link to show the concept:

Update:

Here is a bit of code to show specifically the bit of code I am struggling with. I can currently select a polygon by clicking... but want to select by dropping an icon on the polygon.

$(function() {

$( "#cow" ).draggable({ revert: true});
$( "#map" ).droppable({
  drop: function( event, ui ) {

    **\\How do I select a polygon on drop here?**
    alert( "Dropped!" );
  }
});


texts: Problems while opening .tab file in QGIS

texts:
I'm currently struggling with adding an OS 1:25k (.tab) to my workspace in QGIS.
It comes up with the message

  is not a valid or recognised data source

Any ideas?


texts:
I installed the OpenGeo Suite 4.8 on MacOS 10.13. I can open pgadmin3 from the package, as well as PostGIS. When opening the pgShapeLoader I get the following error:

pgShapeLoader cannot be opened because of a problem.
Check with the developer to make sure pgShapeLoader works with this version of macOS. You may need to reinstall the application. Be sure to install any available updates for the application and macOS.

  Process:               pgShapeLoader-bin [8517] Path:
  /Applications/PostGIS
  Utilities/pgShapeLoader.app/Contents/MacOS/pgShapeLoader Identifier:
  org.opengeo.pgShapeLoader Version:               ??? Code Type:
  X86-64 (Native) Parent Process:        ??? [1] Responsible:
  pgShapeLoader-bin [8517] User ID:               501

  Date/Time:             2018-06-04 13:34:00.559 +0200 OS Version:
  Mac OS X 10.13.4 (17E199) Report Version:        12 Anonymous UUID:
  9D650379-9A83-3E81-917A-0F25D751B9FA

  Sleep/Wake UUID:       1F5D21D1-952B-4A1E-92BC-3CFBF24AC271

  Time Awake Since Boot: 100000 seconds Time Since Wake:       1100
  seconds

  System Integrity Protection: enabled

  Crashed Thread:        0

  Exception Type:        EXC_CRASH (SIGABRT) Exception Codes:
  0x0000000000000000, 0x0000000000000000 Exception Note:
  EXC_CORPSE_NOTIFY

  Termination Reason:    DYLD, [0x4] Symbol missing

  Application Specific Information: dyld: launch, loading dependent
  libraries DYLD_LIBRARY_PATH=/Applications/PostGIS
  Utilities/pgShapeLoader.app/Contents/Resources/lib

  Dyld Error Message:   Symbol not found: _inflateValidate   Referenced
  from:
  /System/Library/Frameworks/ImageIO.framework/Versions/A/Resources/libPng.dylib
  Expected in: /Applications/PostGIS
  Utilities/pgShapeLoader.app/Contents/Resources/lib/libz.1.dylib  in
  /System/Library/Frameworks/ImageIO.framework/Versions/A/Resources/libPng.dylib

Does anybody understand what the error message is trying to tell me? Is there information in it that might help me to open pgShapeLoader?


texts:
I export bunch of DEM files (10000x5000) from Photoscan and I need to combine them with gdal_merge.py but the problem is that it consumes all the memory availible and can't finish processing.

What workaround can be applied to limit RAM usage?


texts: How to select and highlight feature by attribute

texts:
I have the below code which does a select by shape in my ArcMap how do I achieve this with attributes?

1. The map has only one layer
2. There is a filter for the attributes on features

Below is the code:

private void ZoomToDGVSelection(int selectedRightID)
    {
        IMxDocument pMxDoc = ArcMap.Document;
        IMap pMap = (IMap)pMxDoc.ActiveView;

        ESRI.ArcGIS.Carto.ILayer layer = GetLayersClass.GetFieldBoundaryLayer;
        if (layer is ESRI.ArcGIS.Carto.IGroupLayer)
        {
            ESRI.ArcGIS.Carto.IGroupLayer groupLayer = layer as ESRI.ArcGIS.Carto.IGroupLayer;
            ICompositeLayer pCompositeLayer = layer as ICompositeLayer;
            int layers = pCompositeLayer.Count;
            ILayer pLayer = pCompositeLayer.Layer[0];
            IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;
            IFeatureClass pFeatureClass = pFeatureLayer.FeatureClass;
            IQueryFilter pFilter = new QueryFilterClass();
            pFilter.WhereClause = "RightID = " + selectedRightID.ToString();

            IFeatureCursor pFeatureCursor = pFeatureClass.Search(pFilter, false);
            IFeature pFeature = pFeatureCursor.NextFeature();

            if (pFeature == null)
            {
                System.Windows.Forms.MessageBox.Show("This section doesn't exist");
                return;
            }
            pMap.SelectFeature(pFeatureLayer, pFeature);
            IEnvelope pEnv = pFeature.Extent;
            IGeometry pgeom = (IGeometry)pFeature.Shape;
            pMap.SelectByShape(pgeom, null, false);
            pMxDoc.ActiveView.Extent = pEnv;
            pMxDoc.ActiveView.Refresh();
        }
    }


texts: Find points that are unique between two tables

texts:
I have been trying to work out a query that returns the geometry of points that are unique (eg. not supliacte).
Is there a simple way to perform it in just one query? geom is the only common attribute between the tables.
My tables are rather big, about 50000 rows in each.

I need the opposite return to this query:

SELECT point1.geom, point2.geom
FROM schema.point1 AS fkb2012 , schema.point2
WHERE ST_Equals(point1.geom, point2.geom) ;


texts: Combining two layers and removing their overlap using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.0.

I have two polygon layers showing storm surge and I want them both to be combined to look like a smooth transition into each other. However, they have an overlap and it makes the appearance of them poor if they are combined using merge or union.

Would it be easier to convert them into rasters and perform a calculation through that?


texts: How can I know why an OSM way cannot be used for routing, and how can I fix it

texts:
Background

I use osm2po in an application that needs to calculate distances, times and path between locations.

In short, I take data from Geofabrik as input to osm2po to output a directed graph with OSM Nodes as the graph nodes, and OSM Ways as the graph edges.

Later, to calculate distance, time, and path between two locations, I use this graph to locale the nearest nodes, then use one of the osm2po router to build a route between those two nodes.

This is working fine most of the time.

And when a routing can not be done, it is usually because the data from OSM (particularly OSM ways) contains error or is absent.

More to the point

Through debugging a particular case, I located two nodes that cannot be routed from one to the other.

This one and this one.

As you can see, the only way those two nodes can ever be routed together (using a car) is somehow exiting their local streets and going through this way.

After more debugging, I see that osm2po never use this way, I guess it does not see it as a valid segment.

My questions

Looking at the OSM data, how can I know (understand) why this way is not valid for routing ?
Related to my first question: how can I fix (submit modifications to) OSM to allow routing using this way ?

EDIT:

I forgot to say that if I try to route from one node to an other node that are part of this same way, the routing is obviously working fine.


texts: Populate table by randorm values from a list in QGIS

texts:
Assume I have the following list:

 foo = ['a', 'b', 'c', 'd', 'e']

what is the best way to fill a field by random values from above list?
I've written following code, but it doesn't work.
1. Adding this code to function tab:

import random
foo = ['a', 'b', 'c', 'd', 'e']
def randomword():
    return random.choice(foo)

2. Calling above function in Expression section:

randomword()

QGIS Says: Expression is invalid.

Parser Error:
syntax error, unexpected '(', expecting $end

Eval Error:
No root node! Parsing failed?


texts: Exporting layer extent information for 400 layers in QGIS?

texts:
I have over 400 layers, and I need to extract the extent information into a list.

Is there any way I can extract the extent for all 400 layers into a single file?

Copy and pasting one by one isn't practical.


texts: Counting in for-loop using ArcPy Cursor?

texts:
When using the UpdateCursor, is there a "built-in" way to count within the loop?

In the moment I am using a variable i and add 1 after each iteration:

cur = arcpy.da.UpdateCursor("points.shp","ID1")
i = 0
for row in cur:
    row[0] = i
    cur.updateRow(row)
    i = i +1

Is there a smarter way, something like row.count()?


texts:
I am very new at Python,I'm not sure if this is even the appropriate way to do what I want to do, but basically, I'm looking to eliminate all of the small polygons in a feature class that was created to show all of the watershed basins of an area. I will need to run the Eliminate tool multiple times to eliminate all polygons smaller than an acre (since sometimes the tool will combine 2 small polygons and they still end up being less than an acre). Right now I can run the eliminate tool once with:

arcpy.management.MakeFeatureLayer(basinFC,basinLYR)
arcpy.management.SelectLayerByAttribute(basinLYR,"NEW_SELECTION","Acres < 1")
arcpy.management.Eliminate(basinLYR,basin_elim)

Since I've never used a loop in Python before, I'm not sure which to use or where to place it.


texts:
I try to merge 2 raster tiffs in QGIS 3.10 (Raster-Miscellaneuous-Merge). The merge is executed but the resulting map is just black, the scale is set to 0, i.e. all the information of the rasters is gone. I attach screenshots of the configuration of the merge, the rasters before and after the merge and the log, to illustrate my problem.

What mistake have I made?


texts: Problem with repo qgis.org/debian "does not have a Release file"

texts:
When I include

deb https://qgis.org/debian/ buster main
deb-src https:://qgis.org/debian/ buster main

in my /etc/apt/sources.list file, and run

sudo apt uptdate

I get:

E: The repository 'https://qgis.org/debian buster Release' does not have a Release file.
N: Updating from such a repository can't be done securely, and is therefore disabled by default.
N: See apt-secure(8) manpage for repository creation and user configuration details.

This does not happen with

deb https://qgis.org/debian-ltr/ buster main
deb-src https://qgis.org/debian-ltr/ buster main

Is there anything I can do in my side?

My system:

alobo@Delia:~$ uname -a
Linux Delia 4.16.0-2-amd64 #1 SMP Debian 4.16.12-1 (2018-05-27) x86_64 GNU/Linux
alobo@Delia:~$ cat /etc/os-release
PRETTY_NAME="Debian GNU/Linux buster/sid"
NAME="Debian GNU/Linux"
ID=debian
HOME_URL="https://www.debian.org/"
SUPPORT_URL="https://www.debian.org/support"
BUG_REPORT_URL="https://bugs.debian.org/"


texts: How to move projects from a server to local hard disk or removable drive?

texts:
I have a QGIS project on a desktop that I need to move to a laptop to take into the field.

Whilst in the office, both PCs talk to the project on server via a persistent mapped drive (L:) - ie. project lives in L:\My Project\QGIS.

When I go to the field, I mirror the project onto a USB 3.0 thumb drive but in order to access it I need to remap the L drive on the laptop to point at the thumb drive rather than the network server.

This is relatively straightforward but does involve having to remap the L drive back to the network when I arrive back.

I would prefer to tell QGIS that whilst in the field my project is another drive (eg. M) so I can leave the L mapping alone.

In my case all my files live in the same project area so project reparenting would be relatively simple from a coding standpoint. But having said that from QGIS's standpoint, there is no necessity for all files to live on the same drive and thus reparenting would be more complex.

Thinking about what I have just written, perhaps the good old DOS SUBST command is the easiest option for doing what I want to do.


texts: How to find if two gps tracks were from the same trip?

texts:
Given I have gps tracks from many users running or hiking, how can I find out who were running together? In other words, which algorithms are there to find out who were in the same party?

Haven't worked much with GIS so I might use the wrong terms. Matching GPS tracks is almost what I want but it doesn't considering time. My own algorithm would be something like: if two tracks are within 10 meters in distance for 50% of the time they are considered traveling together. Of course the actual tolerances might change. That is easy to code but if there are more general algorithms that calculates the total "distance" or correlation between to gps tracks I would like to use those.


texts:
I have a mesh ( with points and triangles that define it), and I will need to create a DEM model out of it ( in geotiff format if possible for space consideration, if not, then in ASCII grid format). I need the DEM for watershed and catchment delineation purpose.

Instead of using existing software like ArcGIS, I want to write my code to do it ( in C++ or C#), so that I can later integrate them into my own GIS software. is there anything in GDAL or other open source source code libraries that do this well? Or do I have to write from scratch?

PS:

there are a lot of posts like Seeking tool to generate Mesh from DTM? on how to do it the other way round, ie: from DEM grid points to mesh, but they are not what I want.
I have a similar question as Getting grid points for DEM purpose that is bounded by surveyor points? on how to get DEM grid points from point clouds, but that is fundamentally different from this one. I don't need the mesh interpolation part for this question, as I already have the mesh.
I can always create the bounding box for the mesh, loop through each grid points, get the z, and then construct the grid, but this is crude and can be slow. I'm unsure whether there is a faster way of doing it?


texts:
I'm trying to assign weights based on intersection of two layers. Below is my query

SELECT count(*), -- total number of grids in a postcode
       CASE WHEN st_intersects(centroids.geom, urban.geom)
            THEN 0.90::NUMERIC / count(*)
            -- when there is an intersection, 90% weights is equally divided among ALL the intersected grids
            ELSE 0.10 / count(*)
            -- When no intersection, 10% weights is equally divided among non intersected grids
       END  AS Weights
  FROM vrg.gfk_2016_id_5_digit_pcd_areas2013_projected AS postcodes,
       vrg.rdsid86_quadgrid                            AS fishnet,
       vrg.rdsid_86_quadgrid_centroids                 AS centroids,
       vrg.rdsid86_katrisk_poly_projected              AS urban
 WHERE fishnet.gid = centroids.gid
   AND postcodes.id = '66371'
 GROUP BY centroids.geom, urban.geom

The query takes long time to run. How can I modify so that I can create a new table with vrgid and the correspong weights?


texts: Delete single raster from Raster Dataset

texts: Is there any possibility to delete single raster previously loaded to a Raster Dataset?


texts: Why does raster2pgsql change the underlying raster properties?

texts:
I am using raster2pgsql to import a raster of the UK (which involves 2 tiles) from the Copernicus DEM that is in SRID = EPSG 3035.
If I just directly pull these tiles into QGIS, then I find that the raster extremes are:
Min = -2.6m &
Max = 1299.8m
These are reasonable numbers given the highest point in the UK is 1345m and accounting for some error.
When I import these using
raster2pgsql -s 3035 -I -C -M *filepath* -t 200x200 public.raster | psql --DB connection

then I find the extremes are:
Min = -176.4m &
Max = 1524.2m
I've looked to try and understand what exactly raster2pgsql is doing, but I can't find anything on the accuracy. Can anyone help me understand why there is such a difference between the original raster and the imported data in raster2pgsql?
For ref, if I Check the constraints, I get the following but I'm not sure if there are any clues in there or not:
Check constraints:
    "enforce_height_rast" CHECK (st_height(rast) = 200)
    "enforce_max_extent_rast" CHECK (st_envelope(rast) @ '0103000020DB0B000001000000050000000000000060E346410000000060E346410000000060E3464100000000D01253410000000080844E4100000000D01253410000000080844E410000000060E346410000000060E346410000000060E34641'::geometry) NOT VALID
    "enforce_nodata_values_rast" CHECK (_raster_constraint_nodata_values(rast) = '{-340282346638529000000000000000000000000.0000000000}'::numeric[])
    "enforce_num_bands_rast" CHECK (st_numbands(rast) = 1)
    "enforce_out_db_rast" CHECK (_raster_constraint_out_db(rast) = '{f}'::boolean[])
    "enforce_pixel_types_rast" CHECK (_raster_constraint_pixel_types(rast) = '{32BF}'::text[])
    "enforce_same_alignment_rast" CHECK (st_samealignment(rast, '0100000000000000000000394000000000000039C00000000060E346410000000080844E4100000000000000000000000000000000DB0B000001000100'::raster))
    "enforce_scalex_rast" CHECK (round(st_scalex(rast)::numeric, 10) = round(25::numeric, 10))
    "enforce_scaley_rast" CHECK (round(st_scaley(rast)::numeric, 10) = round(- 25::numeric, 10))
    "enforce_srid_rast" CHECK (st_srid(rast) = 3035)
    "enforce_width_rast" CHECK (st_width(rast) = 200)


texts:
My client wants a map of bus routes rendered as white lines with blue casing.  I'm using QGIS 2.8.2 on Mac OSX 10.10.5.  I found a .kml dataset here: http://developer.trimet.org/gis/data/tm_routes.kml.  I also tried using the .shp.

When I add the layer to the map and apply a cased line styling, some of the routes render as nice unified lines and others render as segmented, like earthworms:

I could easily cheat and render a duplicate layer (featuring only the thin white lines) on top of this one.

Why do the lines render in different ways?

Is there a way to render them all as unified?


texts: Using geometrygenerator and different CRS in QGIS?

texts:
In a QGIS3-project (QGIS Version 3.4.3) I'd like to display points from a textfile and constructed lines from this points in an given length and an angle from the textfile with the geometrygenerator.
The CRS of the map is ETRS89 32N (EPSG:25832) and the coordinates from the textfile are WGS84 (EPSG:4326). "On-the-fly-transformation" is turned on.
In the textfile there are the values "x", "y", "xangle".

To display the lines in the given direction, I used this code in the geometrygenerator:

make_line($geometry,
 make_point(
 $x + 0.002*cos(radians(90 - "xangle")),
 $y + 0.002*sin(radians(90 - "xangle"))
 ))

The problem is, that the resulting line differs in the angle from the same construction with points which are original in EPSG:25832 and therefore are in the same projection like the map.
Here one screenshot (blue line = map crs 25832; red line = crs 4326) with 90°,
which looks not to problematic:

Here another one, with angle 152°, which is much more worse:

The other thing is, that ETRS is metric, so I can define the length of the line in Meters and not (like above) in Degrees.
I tested to work in the geometrygenerator with "make_point" direct from the columns of the textfile and "transform". But then there was no line generated.

Maybe someone has got a hint for me to solve this problem. First export the txt to shape with transformation to ETRS is not an option, because the idea is to read from the file in realtime, when something is changed.


texts:
I originally posted this as Install Opengeo 3.0.1 in ubuntu 12.04, but I thought it would be best to ask a new question .  If I should have simply posted on the other forum, my apologies.

I am also having issues installing OG on Ubuntu 12.04.

When I run:

apt-get install opengeo

I get the following:

Reading package lists... Done
Building dependency tree
Reading state information... Done
Some packages could not be installed. This may mean that you have
requested an impossible situation or if you are using the unstable
distribution that some required packages have not yet been created
or been moved out of Incoming.
The following information may help to resolve the situation:

The following packages have unmet dependencies:
opengeo : Depends: opengeo-server (>= 4.0.0) but it is not going to be installed
E: Unable to correct problems, you have held broken packages.

This is my first time using Linux.


texts: Problem when converting sosi to shapefile with ogr2ogr on Linux

texts:
I am trying to convert a SOSI file (Norwegian geodata format) to a shapefile. I have set up gdal, but apparantly there is a lot of things that I dont have support for. When I do a ./configure, this is the result:

  LIBTOOL support:           yes

  LIBZ support:              external
  LIBLZMA support:           no
  GRASS support:             no
  CFITSIO support:           no
  PCRaster support:          internal
  LIBPNG support:            internal
  DDS support:               no
  GTA support:               no
  LIBTIFF support:           internal (BigTIFF=yes)
  LIBGEOTIFF support:        internal
  LIBJPEG support:           internal
  12 bit JPEG:               yes
  12 bit JPEG-in-TIFF:       yes
  LIBGIF support:            internal
  OGDI support:              no
  HDF4 support:              no
  HDF5 support:              no
  NetCDF support:            no
  Kakadu support:            no
  JasPer support:            no
  OpenJPEG support:          no
  ECW support:               no
  MrSID support:             no
  MrSID/MG4 Lidar support:   no
  MSG support:               no
  GRIB support:              yes
  EPSILON support:           no
  WebP support:              no
  cURL support (wms/wcs/...):no
  PostgreSQL support:        no
  MySQL support:             no
  Ingres support:            no
  Xerces-C support:          no
  NAS support:               no
  Expat support:             no
  libxml2 support:           yes
  Google libkml support:     no
  ODBC support:              no
  PGeo support:              no
  FGDB support:              no
  MDB support:               no
  PCIDSK support:            internal
  OCI support:               no
  GEORASTER support:         no
  SDE support:               no
  Rasdaman support:          no
  DODS support:              no
  SQLite support:            no
  PCRE support:              yes
  SpatiaLite support:        no
  DWGdirect support          no
  INFORMIX DataBlade support:no
  GEOS support:              no
  Poppler support:           no
  Podofo support:            no
  OpenCL support:            no
  Armadillo support:         no
  FreeXL support:            no
  SOSI support:              no

  SWIG Bindings:             no

  Statically link PROJ.4:    no
  enable OGR building:       yes
  enable pthread support:    yes
  enable POSIX iconv support:yes
  hide internal symbols:     no

And doing the command results in this:

~/gdal$ ogr2ogr -f "ESRI Shapefile" /media//Geodata/testing/sosi.sos /Documents/test.shp
FAILURE:
Unable to open datasource `/Documents/test.shp' with the following drivers.

I am fairly inexperienced with this, so I am trying to learn. What is going on here, what am I doing wrong?


texts:
I have asked similar question at stackoverflow but posting here as well since it involves ArcGIS and python Add-In. Following code shows first step where i have to declare output folder as global so that later outputs can be saved in it as well. Right now I am getting an error at output folder string r'optfile/ras1'. Any help how to correctly store files in output folder and declare it as global would be appreciative.

import arcpy
import os
import pythonaddins

from datetime import datetime

now = datetime.now()
month = now.month
year = now.year

optfile = "C:/temp/"+str(year)+"_"+str(month)

class DrawRectangle(object):
"""Implementation for rectangle_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.cursor = 1
        self.shape = 'Rectangle'
        os.makedirs(optfile)

def onRectangle(self, rectangle_geometry):
    """Occurs when the rectangle is drawn and the mouse button is released.
    The rectangle is a extent object."""

    extent = rectangle_geometry
    arcpy.Clip_management(r'D:/test',
                          "%f %f %f %f" % (extent.XMin, extent.YMin, extent.XMax, extent.YMax),
                          r'optfile/ras1', "#", "#", "NONE")
    arcpy.RefreshActiveView()


texts: Replacing part of a tileset with gdal2tiles

texts:
Is it possible to run gdal2tiles to replace a specific part of a tileset?
We've created a tileset whose generation took almost 2 days and then found out that one of the input files had an issue (it shows up white). It would be awesome to just replace that one input file and have it tiled so it fits into the final tileset.
I am almost positive I read that gdal2tiles can do that, but I really cannot find anything about it anymore. Was I reading about a different software, perhaps?


texts:
To avoid the possibility of opinions, let me rephrase the question.  Which software packages are able to produce an Oblique Mercator map between any two arbitrary points, similar to what is offered in Matlab (see http://www.mathworks.com/examples/mapping/mw/map-ex21842716-create-long-narrow-oblique-mercator-projection) but without Matlab's price tag of thousands of dollars?

From what I've learned, Google Maps does not seem to work well at the poles.  I'm willing to pay a few hundred dollars for the package if necessary, or I can put in the time to learn to use Geotools, GMT or whatever tool that can produce these types of maps.


texts: Different values while calculating distance

texts:
I'm developing a project using LeafleatJS and PostGIS. Currently I'm getting trouble while calculating distance between two coordinates. For example:

  Point 1 - lat: -12.99835864475412, lng: -38.506194949150085

  Point 2 - lat: -12.999215865191118, lng: -38.50590527057648

Leaflet calculates 99.906 m, while using the following SQL in PostGIS:

SELECT ST_Distance_Sphere(
ST_GeomFromText('POINT(-12.99835864475412 -38.506194949150085)', 4326),
ST_GeomFromText('POINT(-12.999215865191118 -38.50590527057648)', 4326));

I get 81.248 m.

Leafleat uses Haversine distance formula to calculate distance and I'm using its default CRS - EPSG3857 - while on PostGIS I'm storing data as a Geographic Point with SRID 4326.

I know, it seems confusing but I'm following the answer given on this link:

  All of this further confused by that fact that often even though the map is in Web Mercator(EPSG: 3857), the actual coordinates used are in lat-long (EPSG: 4326). This convention is used in many places, such as:

    In Most Mapping API,s You can give the coordinates in Lat-long, and
    the API automatically transforms it to the appropriate Web Mercator
    coordinates.

I'm new on GIS so... why there's a difference between the results (I know they use different formulas to calculate, but almost 20m difference is pretty considerable)? Am I storing and retrieving data in the wrong way? Which one is the most accurate (Leafleat or PostGIS)?


texts: Auto-reprojecting (st_transform) SRIDs in Postgis

texts:
I have a postgis table that is edited from different sources.

The points added from Source 1 are in SRID 23030.
The points added from Soruce 2 are in SRID 4230.

I want all the points in the table to be stored in SRID 23030.

I can create a trigger and st_transform(the_geom,23030) all NEW.the_geom but i was wondering is this is the easiest way to achieve this.

Maybe i can create a CONSTRAINT or something that reprojects (st_transforms) every point if SRID it's not in 23030 or something like that?

Any suggestion?

Regards,


texts: Is there a way to turn on multiple layers at once? I have around 100 shapefiles which are turned off by default on entry. Turning the folder on/off makes no difference and you can’t select multiple layers using Ctl then turn them all on and off.


texts: Selecting not overlapping polygons using QGIS

texts:
I have a shapefile with 1000 points and a layer with 30 000 polygons. Not all of the points overlap with the polygons. For all of these points, which are not overlapping with a polygons, I need to find a polygon, which is closest to the point.
How can I do this spatial join using QGIS?


texts: Using insert and search cursors on same feature class to clone row?

texts:
I'm using the code below in order to attempt to clone a row in a feature class and add that same row to the same feature class after updating its metadata (Name, date_start, date_end, etc.)

The problem is that it is duplicating the rows several times, when expecting 2 rows i get 4, expecting 3, getting 8.

inFeatureClass  = "C:\TS.gdb\IN_TSGIS"
arcpy.MakeFeatureLayer_management(inFeatureClass, temp_layer_name)
insert_feature_table(temp_layer_name, temp_layer_name, query_single_col('UID',UID), update_fields)

def query_single_col(col, value):
    return '"'+str(col)+'" = ' + str(value)

def insert_feature_table(inputfeatureclass, insertfeatureclass, expr, fields):
     with arcpy.da.InsertCursor(insertfeatureclass, ("*")) as i_cursor:
        with arcpy.da.SearchCursor(inputfeatureclass, ("*"), where_clause=expr) as s_cursor:
            for s_row in s_cursor:
                i_cursor.insertRow(s_row)
                break

Any ideas as to why that is?


texts:
I'm trying to normalize a DEM raster with focal and a function to set values between 0 and 1 (iteratingly from the neighboring cells). Thus, for each cell, a 3x3 matrix is to be created and the function applied.

While it works well with the built-in functions (min, max, mean...), I can't get it to work with a normalizing function, ie (x-min(x))/max(x)-min(x)).

Example:

# create a raster
x <- raster(ncol=3,nrow=3)
x[] <- c(13,14,23,46,25,18,32,39,33)
extent(x) <- extent(c(-1,1,-1,1))

#applying the function
res <- focal(x, w=matrix(1,3,3), fun = function(x){(x-min(x))/(max(x)-min(x))})

> res
class      : RasterLayer
dimensions : 3, 3, 9  (nrow, ncol, ncell)
resolution : 0.6666667, 0.6666667  (x, y)
extent     : -1, 1, -1, 1  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
source     : memory
names      : layer
values     : 0, 0  (min, max)

I would expect the central cell to return a value of 0,3636... But it just gives me 0. Any idea what I'm doing wrong?


texts: Why has the total area of polygon been reduced after using the splitting tool in QGIS 3.16.6?

texts:
I am using the splitting tool to split an original polygon (survey boundary) into smaller polygons (vegetation types).
When I sum the total area of the shapefile after splitting (all the veg types) it is less than the area of the original polygon by 0.1 ha (100.61ha vs 100.71 ha).
Is this avoidable as it creates issues when reporting on the data. As I have only split the original polygon into component parts I would assume the total would be the same.
P.S. It has made me wonder if it is an artefact of the boundaries resulting from the splitting. i.e. do the boundaries take up 'space' resulting in less total area in the survey area.


texts: Is there a way to import a geodatabase into a sde on the server and keep the sde domains, while adding the new geodatabse domains?