texts: Point inside st_dwithin with wrong results

texts:
I want to get all points, that are in the same area. Therefore, I have the following query for my table check:

SELECT l.location_id, l.location
FROM check l
INNER JOIN check l2 ON ST_DWithin (l.location, l2.location, 500)
                    AND l.location_id <> l2.location_id;

The column location is a Geometry(Point) data type.

The execution gives me an error and gives me more than 2 million result sets...from initial 1560 rows. The points are normally distributed through Europe and setting the threshold to 500 (m) seem to be fine. Is it may be a prolem with the data type?


texts: Selecting polygons based on geometry (specifically right angles) using ArcGIS Desktop?

texts:
I have some land registry data and I would like to know if there is a way of selecting all the polygons that contain a right angle.

I'm using (basic license) Arcgis 10.5.1.

Here's and image of part of the data with the an example polygon I would like to capture


texts: How to Use ST_MakeLine query in leaflet to connect to Carto

texts:
I want to run in leaflet a query that works in carto Builder to connect points in a line:

SELECT ST_MakeLine(the_geom_webmercator) as the_geom_webmercator, min(cartodb_id) as cartodb_id FROM my_table

i try to reproduce this in an html page :
https://jsfiddle.net/2w4a3f1x/

var sqlQuery = "SELECT ST_MakeLine(the_geom_webmercator) as the_geom_webmercator, min(cartodb_id) as cartodb_id FROM my_table";

// Set CartoDB Username
var cartoDBUserName = "admin";

// Get CartoDB selection as GeoJSON and Add to Map
function showAll(){
  if(map.hasLayer(pointsincartodb)){
    map.removeLayer(pointsincartodb);
  };
  map.setView(new L.LatLng(40.8722012,14.3036309), 10);
  $.getJSON("https://spatial.query.city/user/"+cartoDBUserName+"/api/v2/sql?format=GeoJSON&q="+sqlQuery, function(data) {
    pointsincartodb = L.geoJson(data,{
      onEachFeature: function (feature, layer) {
       layer.bindPopup('<pre>'+JSON.stringify(feature.properties,null,' ').replace(/[\{\}"]/g,'')+'< /pre>');
       layer.cartodb_id=feature.properties.cartodb_id;
      }
    }).addTo(map);

But the link to geojson:

gives an error: {"error":["column \"the_geom\" does not exist"]} but i know that the_geom exist!
as you can see here: https://spatial.query.city/user/admin/dataset/my_table


texts:
I have ArcSDE for SQL Server 2008 R2.

I want to install ArcSDE for another version of SQL Server.

I want to download the ArcSDE software but I cannot find it in search.

Where can I find it?


texts: Extracting polyline segment from points using QGIS?

texts:
I have a polyline from a recent hike going up and down mountain summit.  I would like to extract just the portion of this line that has the track points overlaid.

Is there away to do this without having to trace the portion of this polyline?

Also, I’m trying to do this project strictly in QGIS.


texts: Line style with alternating dots and long dashes

texts:
I´m trying to create that personalized style for hidrography representation (in blue at the Picture); the lines must have 0.2 mm of edge, 3 mm of length, the space between lines must be 1.5 mm and the points, with 0.4 mm of size, must be between in the center of all between-lines spaces. Anybody could teach me how to build that style?

So far, i did this in the Picture. The first point is in the right place. After that everything is in the wrong place:


texts:
My colleague and I are georeferencing a bunch of .tif images.  After each image, we'll click "update georeferencing".  The resultant .tif files are different between our computers and I'm trying to learn why.

If I use my colleague's computer to georeference the file, I can keep the georeferenced .tif and delete the .tfwx, tif.aux.xml, and .tif.ovr files, and the .tif will still appear in the correct location.  ArcGIS recognizes the correct projection.  But if I georeference the file on my computer and delete the .tfwx etc. files, the georeferencing information from my file is completely lost and it does not appear in the correct location.  The coordinate system information is lost on my .tif file.

I'm flabbergasted about why our computers are performing this task differently.  I've followed the exact same steps on the exact same files on both computers.  Can you offer any advice?  We're both using ArcGIS 10.

I've run gdalinfo on the various geotiffs, and the results look identical between the two resultant .tif files (right two panels).  See gdalinfo results here - https://dl.dropbox.com/u/17669584/error%20messages/stav_vs_matt_computer.JPG.


texts: Merging selected features in one DBF table (polygons) using ArcPy?

texts:
I'm absolutely new to Python and programming. How can I merge selected features in one DBF table (polygons) using Python in ArcGIS 10.2 for Desktop?

Here's what I need:

Here's what I've got:

      arcpy.Dissolve_management(S2_merge1, S2_dissolve5, "Name5", "", MULTI_PART", "DISSOLVE_LINES")

      arcpy.JoinField_management(S2_dissolve5, "Name5", S2_merge1, "Name5", "Name1")

This works but I need to do this without dissolve command - You see I need to merge selected features in one shp/dbf like I do this in Editor session. I can't find this type of toolbox/command in ArcGIS

And is there any possibility to configure the Dissolve tool to get shp-theme with all fields what was in original shapefile like it was in ArcView 3.2?


texts: Why do we use the terms “longitude” and “latitude” for a spherical body?

texts:
To specify a geographic coordinate we use the terms “longitude” and “latitude”. Far as I know, these terms derive from the Latin latitudo (=
breadth) and longitudo (= length). In German we use the terms “Breite” (=english “breadth”) and “Länge“ (= english “length”).

The earth is a sphere. Why do we use “breadth” and “length” here? The earth has the same breadth and length in all directions. We may also use the terms “longitude” and “latitude” just the opposite.


texts: How to select vector features with partial (>33%) overlap with another vector feature?

texts: I have a grid and a range map, and would like to select those grid squares which overlap at least 33% with the range map.  The range map is a single vector feature.  (I'm new to QGIS but have some Arc experience) Thanks!


texts: Using the ArcGIS Time Slider on non-temporal data?

texts:
We're looking at re-purposing ArcGIS temporal tools, chiefly the Time Slider, to model potential flood hazards. Our input data are a series of polygons showing varying water depth. Elevation (depth) is stored as metres in an attribute field as well as PolygonZ. The field value range is small, e.g. 316.0 to 320.8.

On of the immediate issues we encounter is that the tools require data to be in YYYY/YYYYMM/.../YYYYMMDDhhmmss format and truncate decimals. It would be nice to be able to leverage the automatic labelling too, but that's not essential.

So, how would you restructure the input data to make maximal use of the temporal tools? The primary result will be used in ArcGIS Desktop, but there is talk of publishing on ArcGIS Online as well. So try to accommodate both avenues if feasible.


texts: Create LineString from lat/lon coords form CSV. with beginning and ending lat/lon?

texts:
I am iterating through a CSV file of tornadoes, and for each row getting the beginning lat, beginning long, ending lat, and ending long. and trying to plot their line geometry on a map. I got the points for the starting points but I am having troubles creating the lines to show the path of the tornadoes.

I have:

for i, row in dataframe.iterrows():
    lx = row[begin_long],[ending_long]
    ly = row[begin_lat],[ending_lat]
    line_geometry.append(LineString(lx,ly))


texts: Extraction of Cloud Coverage percentages of Sentinel 2 images

texts:
I have to download S2A images over a large area from 2016 March to 2019 January. Before start downloading data, first I want to know cloud coverage percentages of all the images over this period.

Is there any way to do it either on QGIS or using python code?


texts:
I'm using R with ggplot2 to generate maps from shapefiles and from data in SQL Server. I'm trying to get to a stage where I can automate a large chunk of the mapping I do. So far this is working well and I'm able to generate some fairly nice looking maps (just need to work on labelling).

I'm wondering if it's possible to create a geographical heatmap from a list of point co-ordinates (which are irregularly spaced)? Similar to using geom_density2d, but raster.


texts: Spatial Join Hangs on Two Specific Shapefiles in QGIS 2.8

texts:
I have successfully done spatial joins with QGIS many times. However, with two shapefiles QGIS hangs or never seems to complete. (I’ve let the process run for over two hours.)

I’ve tried all the obvious things that I can find by searching and reading other posts, e.g., making sure the two shapefiles have the same CRS (link), but it always hangs. What’s more confusing is that I can spatially join both shapefiles with the same other shapefile, but I can’t join them directly. IOW, spatial-join(A,B) and spatial-join(B-C) both work, but spatial-join(A,C) does never completes.

Do you have any tips, suggestions, insights about what might be going on?

DETAILS

I have three layers for a state with successively smaller shapes:

Congressional Districts (CDs) — order a dozen shapes
Precincts — order 2,500 shapes
Blocks — order 300,000 shapes

The CRS for all three layers is EPSG:4269, NAD83. However, the CRS for the Precincts and CDs shapefiles that I got from the state was something else. So, I changed the CRS to EPSG:4269, NAD83 using “Save As …” (per this).

I want to assign Precincts to Blocks and CDs to Blocks, using Vector > Data Management Tools > Join Attributes By Location … Again, I’ve done this kind of operation with many shapefiles successfully, so I believe there’s something subtle going on somewhere that I’m missing.

Vector > Data Management Tools > Join Attributes By Location ... works for

Blocks and Precincts layer with “Target vector layer” = Blocks and
“Join vector layer” = Precincts; and
Precincts and CD layer with
“Target vector layer” = Precincts and “Join vector layer” = CDs

However, it hangs, never completes, for

Blocks and CDs layer with “Target vector layer” = Blocks and “Join vector layer” = CDs

The first two spatial joins complete relatively quickly, taking no more than a minute or two. The latter, transitive spatial join never completes though. As I said, I’ve let it run for two hours to no avail.

NOTE: The only reason I spatially joined Precincts and CDs was to see if it would work. (It does.) A fraction of Precincts actually span CDs which is why need to do this transitive join.

In the process of making sure everything was set up right, I made the shapes in the coarsest layer (CDs) transparent, so I could see the smaller shapes (blocks) through them, and the align correctly. So, I believe the coordinate systems are set up correctly.

I’m running QGIS 2.8 on Mac OS X 10.11.3 (El Capitan).


texts: How to aggregate point measurements in areas?

texts:
I have a shapefile with the geographical boundaries of towns in a specific region;
then i have a csv file (with xy coords) with a lot of measurements about pollution in this region.

I'd like to have the mean value of pollution for town, obviously averaging over the measurements inside each town.
How can i do it? I'm new to QGIS.


texts:
Using code from http://resources.arcgis.com/en/help/main/10.1/index.html#//00170000016p000000 trying to enable editor tracking. It works fine if there is no feature dataset in the GDB, but the moment I add a feature dataset in the GDB I get the following error:

[u'QGC_LINK', u'UDM']
Test
Enabling editor tracking on QGC_LINK
Line 60
Failed to execute. Parameters are not valid.
ERROR 000110: QGC_LINK does not exist
Failed to execute (EnableEditorTracking).

It seems that it can't find the file path.

The file structure is as follows:

Test.gdb
Test.gdb/UDM
Test.gdb/QGC_Link
Test.gdb/FD/UDM_1

Code as follows:

import arcpy, os

# Set the workspace
workspace = arcpy.GetParameterAsText(0)

# Set the workspace environment
arcpy.env.workspace = "M:/GIS/Test.gdb/"

# Get all the stand alone tables and feature classes
dataList = arcpy.ListTables() + arcpy.ListFeatureClasses()

print dataList

# For feature datasets get all of the featureclasses
# from the list and add them to the master list
for dataset in arcpy.ListDatasets("", "Feature"):
arcpy.env.workspace = os.path.join(workspace,dataset)
dataList += arcpy.ListFeatureClasses()

print dataset

# Execute enable editor tracking
for dataset in dataList:
print 'Enabling tracking on ' + dataset
arcpy.EnableEditorTracking_management(dataset, "Created", "CreatedDate", "Modified", "ModifiedDate", "ADD_FIELDS", "UTC")
print 'Enabling complete'


texts:
I would like to export with ogr2ogr the data I have in a kml as ExtendedData

        <Placemark>
        <ExtendedData>
            <SchemaData schemaUrl="#15540_punti">
                <SimpleData name="myid">9</SimpleData>
                <SimpleData name="IDENTITA">15540_IX</SimpleData>
                <SimpleData name="NOME">Fontana</SimpleData>
                <SimpleData name="AMP_CONO_V">80</SimpleData>
                <SimpleData name="PROF_1SSP">5,2</SimpleData>
                <SimpleData name="PROF_1P">26</SimpleData>
                <SimpleData name="PROF_2P">78</SimpleData>
                <SimpleData name="PROF_3P">156</SimpleData>
                <SimpleData name="PROF_SFOND">245.26</SimpleData>
                <SimpleData name="ASSE_DIR">117g39p8s</SimpleData>
                <SimpleData name="UNITA PAES">Piazza Santa Croce</SimpleData>
                <SimpleData name="QUOTA_SLM">276 + 1.50</SimpleData>
            </SchemaData>
        </ExtendedData>
        <Point>
            <coordinates>11.26078815736627,43.76921971994336,0</coordinates>
        </Point>
    </Placemark>
    <Placemark>

I have tried with gdal 1.9.2 and libkml to export in CSV, but I obtain a file that does not contain the fields of the original kml (myid, IDENTITA, NOME, etc.). This is the output

Name,description,timestamp,begin,end,altitudeMode,tessellate,extrude,visibility
,,,,,,-1,-1,-1
,,,,,,-1,-1,-1
,,,,,,-1,-1,-1
,,,,,,-1,-1,-1
,,,,,,-1,-1,-1
,,,,,,-1,-1,-1

Is it possible with ogr2ogr to export  KML ExtendedData in tabular form?

Thank you,

Andrea


texts:
I want to get population counts for places in California. [Ideally Census Tracts, but that's not a supported geography in the population dataset?]

After reading through the documentation and requesting an API Key I built this URL:

    https://api.census.gov/data/2016/acs/acs1?get=B01001_001E,NAME&for=place:*&in=state:06&key=<api_key_goes_here>

I get the expected response in JSON format when I type it into a web browser but when I go connect the dataset in Carto, I get an error about supported files. Is the URL formatted properly? Am I connecting it incorrectly? I've been using the "connect dataset" and pasting the URL there.


texts: How to properly establish a PostgreSQL connection using QgsCredentials

texts:
I'm trying to establish a Postgres connection using QgsCredentials as described here:

How to access user credentials in database connection?

The only problem is that I cannot establish a working connection unless I "physically" add the layer to the map canvas (which I don't want because I only want a connection to the DB and a reference to the layer which will be edited later on in the code)

When clicking a button the run function reads:

uri = QgsDataSourceURI()
connInfo = uri.connectionInfo()

(success, user, passwd) = QgsCredentials.instance().get(connInfo, None, None)

if success:
    # Isn't this enough to establish the connection and get a valid reference to the layer?
    uri.setConnection("HOST", "PORT", "DB_NAME", user, passwd)
    uri.setDataSource("SCHEMA", "TABLE_NAME", "the_geom")
    LYR = QgsVectorLayer(uri.uri(), "LYR", "postgres")

    # Unless I add the following line, I am not able to use properly the object
    # (which asks the user the credentials until he/she inserts the right ones and gets
    # a working connection with a valid reference to the layer)
    QgsMapLayerRegistry.instance().addMapLayer(LYR)

What am I doing wrong?

Thanks in advance,

m.


texts: Struggling with ArcGIS REST API and tokens

texts:
I'm doing my best to secure my AGS (9.3.1) services, but have struggled with tokens. As you can see, my nw account belongs to the agsadmin and agsusers roles:

I have assigned these same roles to my service and its containing folder:

Generating a token is no problem:

But alas, upon logging in to the REST interface, I am presented with this Unauthorized access message:

I found a post on the ESRI forums describing the same problem, but responses are inconclusive. ("Please make sure that the Security settings for this map service is correct.", blah blah.)

Has anyone encountered this? What could I be doing wrong?

Update: we have abandoned AGS-based security because of its many flaws and capricious inconsistencies. Beware!


texts: Connecting Oracle Linux database to AcGIS 10 on windows PC

texts:
I am wondering if anyone has tried to connect from ArcGIS 10, to an oracle database stored on a Linux server?
I have no problems connecting to oracle from Arc when the data is stored on a windows box, but no matter the different combinations of host name / server name I insert to try and connect to the database which is stored on a Linux box I get an error.

Even if its a case of this is not possible, it would be nice to know so I can try and proceed down a different route.


texts: Splitting a RASTER data (a DEM) using Quadtree or R* tree

texts:
I want to create a tree of raster images, of which root node would be the largest which is original\input raster. The child nodes of root node would be the resultant rasters produced after splitting using splitRaster().

The below code splits original raster into 4 children rasters. The problem is I am not able to save these children rasters as nodes of a tree(a quadtree\ R*).

Below is the Rscript which successfully splits raster recursively.

  require(raster)
  require(rgdal)
  require(foreach)
  require(SpaDES)

  b <- raster(file.choose(), package = "raster")
  ra <- b[[1]] # use first layer only
  parent<-0
  xlist<-c()
  splita<- function (r,level,parent){
    tmpdir <- file.path(tempdir(), paste("splitRaster-AtLevel",level))
    strt<-Sys.time()
    nx <-2
    ny <-2
    outpath <- paste("D:/project/outputrasters/AtLevel",level)
    dir.create(outpath)

    # parallel cropping
    if (interactive()) {
      n <- pmin(parallel::detectCores(), 4) # use up to 4 cores
      beginCluster(n)
      y3 <- splitRaster(r, nx, ny, c(0.5, 0.5), path = file.path(outpath, paste("Level",level)))
      files <- list.files(path=tmpdir,pattern = "rasterlayer_tile.*.gri$")
      outfiles <- paste0(outpath, files)
      extension(outfiles) <- 'tif'
      iters=4
      print (paste("splitting raster into 4 regions at level ",level ))

      for(i in 1:iters){
        print(paste("max elevation of region ",i,"is",maxValue(y3[[i]])))
        rc <- writeRaster(y3[[i]],file = paste0(outpath, paste("rasterlayer_tile", i, sep = ".")),bylayer =TRUE ,format="GTiff",overwrite= TRUE)
      }
      endCluster()
    }
    print((level)
    if (parent<4) {
      iters<-4
      for(i in 1:iters){
        splita(y3[[i]],level,parent)
        parent<-parent + 1
      }
    }
  }
  splita(ra,1)
  unlink(tmpdir, recursive = TRUE)

  print(Sys.time()-strt)


texts: Eliminating raster gaps using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.5.

What would be best solution to fill gaps and empty data when forming mosaic?

I tried mosaic, mosaic to a new rater. Not satisfactory.

Is there a way to blend with basemap?


texts: Labeling each polygon based on the number of points inside it using Python

texts:
I have a parcel layer and a trees layer, I want to label each parcel with the number of trees it contains and the sum of DBH. I wrote a script that calculates this and it works stand-alone, but I can't figure out how to get it to work with a Label Expression.
Here is the code I've tried (ArcMap 10.2):

import arcpy
import numpy
def FindLabel():
  mxd = arcpy.mapping.MapDocument("CURRENT")
  parcels = arcpy.mapping.ListLayers(mxd, "Parcel*")[0]
  trees = arcpy.mapping.ListLayers(mxd, "Surveyed Trees*")[0]
  arcpy.SelectLayerByLocation_management(parcels, "CONTAINS", trees, "#", "NEW_SELECTION")
  arcpy.SelectLayerByLocation_management(trees, "WITHIN", parcels, "#", "NEW_SELECTION")
  with arcpy.da.SearchCursor(parcels, ["OBJECTID"]) as cur:
    for row in cur:
      arcpy.SelectLayerByLocation_management(trees, "WITHIN", parcels, "#", "NEW_SELECTION")
      treeCount = int(arcpy.GetCount_management(trees).getOutput(0))
      field = arcpy.da.TableToNumPyArray(trees, ("DBH"), skip_nulls = True)
      sum = field["DBH"].sum()
      message = "Trees = " + str(treeCount) + "\n" + "DBH = " + str(sum)
      return message

This sends ArcMap into a long loop and then seems to return the last calculated value for each parcel, so everything is labeled "1 Tree 7 DBH".

Maybe I shouldn't be using a SearchCursor?

Solution

Here's the code I ended up using. It's still processor intensive, but this fits our workflow fine.

import arcpy
def FindLabel ( [Loc_Name] ):
  mxd = arcpy.mapping.MapDocument("CURRENT")
  trees = arcpy.mapping.ListLayers(mxd, "ProgramSurveyGround*")[0]
  arcpy.MakeTableView_management(trees, "treesTable")
  parcels = arcpy.mapping.ListLayers(mxd, "ParcelProcessing*")[0]
  Query = ' "Loc_Name"' + " = '" + [Loc_Name] + "'"
  treeSelect = arcpy.SelectLayerByAttribute_management("treesTable", "NEW_SELECTION", Query)
  treeCount = int(arcpy.GetCount_management("treesTable").getOutput(0))
  dbhArray = arcpy.da.TableToNumPyArray("treesTable", ('DBH'), skip_nulls = True)
  dbhCount = dbhArray['DBH'].sum()
  myText = "Trees : " + str(treeCount) + "\n" + "DBH : " + str(dbhCount)
  return myText


texts:
Let's say I have two WFS layers from different sources. Both layers have the same data structure and application schema.

Is there an open source software package like Geoserver that I can use to "merge" them into one WFS layer which other services can then consume?

I've looked at cascading WFS but it doesn't seem to be what I need.


texts: MapBox GL and styling for data driven labels

texts:
I can use the data driven feature to color points, works great

MAP.addLayer({
"id": "schools",
"type": "circle",
"source": "markers",
"layout": { },
"paint": {
    "circle-radius": 8
    'circle-color': {               // for data driven colors need v0.18.0
        'property': 'type',
        'type': 'categorical',
        'stops': [
            ['ELEMENTARY', '#009900'],
            ['MIDDLE', '#cc00ff'],
            ['HIGH', '#1666ff']
   }
}
});

I also want the label to be the same color as the point. All my text labels default to the first stop color. Can I use the data driven feature for text labels, and if so what am I doing wrong?

MAP.addLayer({
"id": "school_label",
"type": "symbol",
"source": "markers",
"layout": {
    "text-font": ["Open Sans Regular"],
    "text-field": "{full_name}",
    "symbol-placement": "point",
    "text-size": 10
},
"paint": {
    //"text-color": "#ff0000",
    'text-color': {
        'property': 'type',
        'type': 'categorical',
        'stops': [
            ['ELEMENTARY', '#009900'],
            ['MIDDLE', '#cc00ff'],
            ['HIGH', '#1666ff']
    }
    }
});


texts: How can I load and draw predefined chart images (diagrams) rule based and related to a point vector layer with QGIS

texts:
I have a bunch of diagrams for a set of route segments, stored as PNG, PDF or SVG with a regular name (topic-segment-id.(png|pdf|svg) for example), representing a certain aspect (like crossing rate or and angle for ships at a segment (..see graphics, polar density plots).

I'm able to generate the images of the diagrams (including a subsequent image processing chain) with a name carrying the segment-ID. I can also create a point layer for the placement of the diagrams, with  attributes for the filename and/or the segment-ID,  at the map (..see map points 1 to 17 corresponding to the diagrams above)  along the route.

How can I load and place the segment related images (the generated diagrams for example) at the addressed points/segments for the corresponding map under QGIS wit the print composer/ map layout tool. I want to have something like that:

PS: I use the statistic language R to create diagrams and  the point layer (package sf) and QGIS to generate the maps.


texts: Debugging code to export .shp that gives invalid SQL statement was used?

texts:
I have a shapefile composed of a number of polygons. The shapefile has table column entitled as name (field: NAME).

How can I export with ArcPy each polygon so I get a shapefile for each one and it's name (ex: name.shp) corresponds with the ones in the table column?

I've attempted to use the following code:

import arcpy
shp = "c:/temp/folder/shapefile.shp"
TargetFolder = "c:/temp/folder"
query = """"NAME" > 2"""
with arcpy.da.SearchCursor(shp, ["NAME","SHAPE@"], query) as cursor:
    for row in cursor:
        arcpy.FeatureClassToFeatureClass_conversion(row[0], TargetFolder, row[1])

I get the error code: An invalid SQL statement was used.

I used code from Export to SHP with ArcPy with query as the basis for mine.


texts: Performing Vector Affine Transformation using QGIS?

texts:
Due to an incorrect datum coordinates, I need to move a large number of elements in shapefiles by a specific distance. In 2.8 I would have used the Affine Transformation tool, but this does not appear to be available to QGIS 3.6.

Is there another tool or process I can use?


texts: Pixel Size field in gdalinfo command meaning

texts:
What is the meaning of the Pixel Size field in the gdalinfo output on a TIFF image?

I have 2 images, one is a PAN image (size 400MB) while one is a DEM (size 25MB). Both are in TIFF format. When I run the gdalinfo command on these 2, I get the following:

For the PAN image:

Pixel Size = (0.600000023841858,-0.600000023841858)

For the DEM:

Pixel Size = (0.000277777777778,-0.000277777777778)

I also get the same values (with one extra value) when I run another command, exiftool, on these 2 images:

For the PAN image:

Pixel Scale                     : 0.600000023841858 0.600000023841858 0

For the DEM:

Pixel Scale                     : 0.000277777777778 0.000277777777778 0

I can see that the output is same, except that gdalinfo puts a minus sign before second value, and exiftool gives 3 values where the 3rd value is 0.

What should I infer from these values? Which image has a better spatial resolution? I know its the PAN image that has a better resolution, but how do I understand that from these values?

I say this because since 0.00027 < 0.6, so one might think that the DEM has a better resolution. Just like when we say that an image has a spatial resolution of 0.25 meter while other has 0.3 meter resolution, it means the 0.25 meter one (the lower one) has better resolution. Or maybe I am understanding it wrong if I'm looking at the value of the number in this case.

Whats the unit of these numbers? Is it meters? If not, then what? Which gdalinfo field could tell me the unit? Right now I get this metadata field in gdalinfo which has the word 'unit' in it:

TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)

How to I explain someone that in this context, bigger number means better resolution, not the other way?

Some additional info that I feel might be helpful:

The metadata field AREA_OR_POINT is as follows:

For the PAN: - POINT

For the DEM: - AREA


texts:
How does multiuser access work in QGis 3.0?

Can multiple people open and edit the same object in a layer and save it at the same time? (Edit: Yes, but it may cause loss of data)
Does it warn you before overwriting or does it update automatically after someone edits it? (Edit: No)
I am using shapefile data at the moment.
Does it work with that or do I have to set up a database-server?
If yes, how can I do so?

Main goals of the database:

all users have read/wirite access
data edited by multiple users at the same time gives the second user a notification/warning before saving


texts: Zonal Statistics or Tabulate Area by Multiple, Cascading Zones

texts:
A basic GIS task that I regularly come across in my work is to calculate the area of a variable of interest by multiple, overlapping and/or hierarchical zones. I am looking for a more efficient way to perform this and export to tabular format for inclusion in reports. Can anyone recommend an approach (or ArcPy script) that can calculate land cover type (raster) by two or more overlapping zones, in an efficient manner, e.g.:

Calculate area by Zone 1
Calculate area by Zone 2, including specifying areas that do and do not overlap with Zone 1
Retain fields in attribute table for both zones
Optional: Export to CSV or Excel, or compile multiple tables into a single table if above was iterated. Example below, where management choice is variable of interest to be reported back by zone 1 and zone 2.

This is similar to using Tabulate Area with the administrative area by zone (Tabulate Area Using Multiple Shapefiles and Multiple Rasters), but I need one or more additional levels of zones, and my zones are not in the same input file.

What I currently do (works, but seems less efficient):

Clip and Iterate:
(1) Clip zone 2 by zone 1 (land use clipped by administrative zone, in example above), (2) iterate through each category in variable of interest (e.g. management choice) and run zonal statistics (cell count) or tabulate area, and (3) export table for each to Excel. This works, but reports total theoretical areas, not accounting for the fact that once a management choice is made it precludes other choices if it overlaps in that area.

Combine:
Another approach I have used, but is also not efficient is to use ArcMap’s ‘Combine’ tool, join back text descriptions to the output so I know what combination code (similar to overlap of Management Choice Codes above, e.g. 1101 (wetland, managed wetland, and forest restoration overlap)), and run zonal statistics on this. Then calculate area based on cell size, and export to CSV/Excel. This works for 1-2 zones, but gets complicated when there are multiple potential values for any given location. For example, interpreting what values should join to 1101 can take time.

Other:
I have not tried this approach as I do not have the ability to install QGIS at work. However, SAGA’s Cross-classification and Tabulation tool might be able to do this too, as referenced here (Is there a "Tabulate Area" function in QGIS?) and here. For example, use Zone 2 as ‘Input Grid 1’ and Management Choice as ‘Input Grid 2) in SAGA, then analyze by Administrative Zone (zone 1, above).

Thoughts on any of the above, or can you recommend another approach that would perform this in a single step (or iterated steps) and report the total area of variable (Management Choice), by Zone 2 (Land Use) within Zone 1, by Zone 1, as well as the total overlapping and non-overlapping?


texts: Loading different shapefile in SharpMap?

texts:
I am currently working on sharpmap project with the need to work on offline maps. As i am fresher in this field,I am following the sharpmap tutorial and facing a problem with loading new shape files in the given tutorial code.

For Example :

SharpMap.Layers.VectorLayer("States");
        vlay.DataSource = new SharpMap.Data.Providers.ShapeFile("path_to_data\\states_ugl.shp", true);

At this line of code,if i pass a different shapefile,code builds with a blank background or no display.
I have tried with different shape files with different sizes but the result is the same. It only works for the mentioned states_ugl.shp file given in the code.


texts:
I need to classify tree heights in a DSM of a slope-y terrain. It´s not a forest but several individual / detached trees. In a plane terrain it won´t be an issue but in that uneven terrain I need somehow to extract the height of the single trees as related to the surrounding ground.

Is there any way to do so using QGIS?

The DSM file is a Tiff, but I can also get it as a point model if this helps.


texts: Wrong northing value: standardizing with a correction by 10 000 000 meters

texts:
I'm trying to draw points on a map and I have a problem using the projections. I have no problem getting points from longlat to UTM. But it seems that there is a problem when specifying UTM zone 15M and then try to plot them. I figure out that there is a difference of 10 000 000 meters between the 2 coordinates that I had. Why is that?

library(sp)
library(mapview)

gps.pt = data.frame(lon = -90.22139,
                    lat = -0.6893730)
lonlat= sp::SpatialPoints(gps.pt, proj4string=CRS("+proj=longlat +datum=WGS84"))
spgeo <- sp::spTransform(lonlat, CRS("+proj=utm +zone=15M +datum=WGS84"))
# OK!
mapview::mapView(spgeo)

# Wrong map (the easting and northing values were collected from [memento][1] on an Android smartphone)
gps.pt2 = data.frame(easting = 809289.5,
                     northing = 9923713)
spgeo2 = sp::SpatialPoints(gps.pt2, proj4string=CRS("+proj=utm +zone=15M +datum=WGS84"))
mapview::mapView(spgeo2)

# Now OK!
gps.pt3 = data.frame(easting = 809289.5,
                     northing = 9923713-10000000)
spgeo3 = sp::SpatialPoints(gps.pt3, proj4string=CRS("+proj=utm +zone=15M +datum=WGS84"))
mapview::mapView(spgeo3)

Do I actually have to subtract 10 000 000 meters or I've just specified my CRS wrong? When should I apply a 10 000 000 meters correction?


texts: Iterate st_intersection (sf) over factor levels

texts:
I wrote a code using sf objects in order to count the number of points per category within cells. The code works well, but I would like to iterate over factor levels.
Here is a reproducible example :
library(sf)
library(raster)
library(dplyr)

region <- st_read(system.file("external/lux.shp", package="raster"), quiet = TRUE)
PT <- st_sf(st_sample(region, size =100, type='random'))
PT$Class <- sample(LETTERS[1:4])
Mesh <- st_sf(st_make_grid(PT, cellsize = 0.2))
Mesh$ID <- seq.int(nrow(Mesh))

intersection <- st_intersection(Mesh, PT[PT$Class == "A",])
PT_result_A <- intersection %>%
  group_by(ID) %>%
  count()

The above works well to count the number of points with the "A" label inside each cell. In my original code/data I have this block as many time as I have categories, and I have a lot of categories.
Is it possible to automate the above to obtain the count of "A", "B", "C" and "D" without repeating the last block?


texts:
I want to map a list of facilities around the UK by their Postcodes. I have no Eastings or Northings just postcodes. Is this possible via a Plugin or do I need coordinates.
I am using QGIS Version 3.10.0


texts: Adding custom columns to the details view in ArcCatalog

texts: When you right click an MXD in ArcCatalog and go to properties you have the option to add a description to your MXD.  The description field can also be set from ArcMap.  What I would like to know is if there is any way to add this description field as a column within ArcCatalog detail view.  There are many ways to customize which columns appear in the detail view using the Customize/ArcCatalog Options/Contents tab.  I've read a few instructions for 9.3 that led me to believe you could do this through the use of XML tags but I haven't been able to achieve anything like this in 10.3.  I understand that MXD's are the only filetype that would contain the property I want in the "description column" and I know the details columns are typically based on properties that would exist in any file type.  Maybe this is a pipe dream, but if anyone has done this please let me know.


texts: Need to quickly know (CLI) whether point lies inside polygon

texts:
in my app I am parsing an XML from the web which gives me such coordinate info:

I am able to use my app to get a pair of float values like 14.871381 and 50.065968. I do not want any complicate processing (using simple C++ program running in CLI) - if possible the best would be if some free library or CLI utility could say whether this point (according to the coordinates) lies in a polygon (in my case a particular region - just one) - maybe like this:

Have a map with defined polygon (region/county) A
Query: 14.871381 ; 50.065968 Reply: A
or:
Query: 14.871381 ; 50.065968 is in region A? Reply: True (or False or 1/0 or something simple)

Most important is to have really efficient and low demanding solution, no GUI needed (program runs in CLI) and it runs on weak HW (single-core Atom CPU) running Kubuntu linux and about 1-2 queries per second.

PS: my brother is able to help me with GIS related information but not with the programming...


texts: Modifying style (line width/color) of vector layer in PyQGIS?

texts:
I'm currently updating some python code from QGIS 2.18 to QGIS 3.4

in 2.18 i used the following code :

#First off I created a virtual layer using an sql query here

vlayer = QgsVectorLayer( "?query=...., "virtual" )

# Then I assign the properties of the polyline layer here

props = { 'width' : '3', 'color' : '255,0,0' }
sl = QgsSymbolLayerV2Registry.instance().symbolLayerMetadata("SimpleLine").createSymbolLayer(props)
s = QgsLineSymbolV2([sl]) vlayer.setRendererV2( QgsSingleSymbolRendererV2( s ) )

I know there have been a number of changes in QGIS 3 for exempleQgsSymbolLineSymbolV has become QgsLineSymboland setRendererV2 which becomes setRenderer. So, with the help of this blog this GitHubthe API and the most recent doc on the QGIS site, I've put the script below together...

layer = QgsProject.instance().mapLayersByName("t_alpha")[0]

single_symbol_renderer = layer.renderer()

symbol = single_symbol_renderer.symbol()

symbol.setColor(QColor.fromRgb(255, 0, 0))
symbol.setWidth(3)

layer.triggerRepaint()

qgis.utils.iface.layerTreeView().refreshLayerSymbology(layer.id())

But this returns an error:

  AttributeError: 'QgsCategorizedSymbolRenderer' object has no attribute
  'symbol'

I've dug through the doc and otherwise tried to remodulate the code a bit but I'm a bit stuck (i'm still learning Python).


texts: How are link costs calculated from osm.pbf?

texts: I imported osm.pbf data into PostGIS using osm2po. The table I obtained (at_2po_4pgr) includes cost data. I would like to know how these costs have been calculated?


texts: ArcGIS Enterprise Geodatabase - MSSQL Server Permisions & Schemes

texts:
I have an active directory group called GIS_Admins. I have created a Server Login that points to this group and the group is mapped to a user for the appropriate database. The User has been made a member of db_owner and the default scheme is set to dbo.

For the most part this setup is working fine. Here is the issue. When a user connects to the Geodatabase using their windows authentication, if the user creates a new feature class it creates a scheme after their domain name. For example the feature class would be named "database.'domain\username'.FeatureclassName". I would like all feature classes created by members of this group to be created into the dbo scheme. What am I missing here?


texts: Latitude and longitude of NetCDF layers in ArcGIS

texts:
I've been having a lot of trouble getting the Multidimension system toolbox to work in ArcGIS (10.3). In particular, the Make NetCDF Raster Layer tool does not pick up on numerical Latitude and Longitude fields, meaning that the tool cannot correctly display my raster.

The sample data I am testing with can be downloaded here, although I've had the same issue with other NetCDF datasets. (There are plenty listed here.)

I have no problems loading this or other datasets with NASA's Panoply, suggesting the issue is entirely with the tool I'm trying to execute, not the dataset. Below, you can see the variables listed, which include longitude and latitude, with type 2D. How can latitude and longitude values have two dimensions on their own?

This still works, as it has worked for every .nc dataset I've tried in Panoply.

However in the Make NetCDF Raster Layer tool in ArcGIS, longitude and latitude cannot be selected
for the X Dimension or Y Dimension. Inappropriately, one (and only one) of these fields can be used as the Variable value. I couldn't get a screenshot with the drop-down showing, but latitude and longitude cannot be specified for the dimension fields. In practice this means that although other values can be chosen for X and Y, none of them give a true geographic representation. That is, the tool cannot actually show me the spatial representation of the NetCDF dataset.

Is there an issue with the Make NetCDF Raster Layer in ArcGIS 10.3? An issue with the dataset? An issue with my workflow? How can I display this dataset in ArcGIS?

EDIT: I have been successful using this tool in ArcGIS with a different NetCDF dataset, for a dataset that Panoply identified the latitude and longitude values as 1D type, rather than 2D type. So I suspect that is the problem. However I am still at a loss as to:

How Panoply seems to get by with the dataset while ArcGIS does not; and
How I can convert this and other datasets where the type of latitude and longitude is (apparently incorrectly) given as two-dimensional, so that ArcGIS can use them.


texts:
I have a problem with Extend Line tool in ArcPy, it cannot extend some lines, i tried to get around by erasing the edges of lines that i want to extend without failure, but despite this some fewer lines remains not extended to any feature in the same layer.

an example of Extend commande

arcpy.ExtendLine_edit("Line_fc" , "", "FEATURE")


texts:
I'm new to the QGIS Pyculator and need to reclass a field ('AP') based on another field ('Mag'). What I have so far, from the ArcGIS field calculator:

def Reclass(Mag,AP):
    if Mag>= 2000:
        return "M"
    elif Mag>80 and Vector_Mag<2000:
        return "L"
    elif Mag<=80 :
        return "S"
    else:
        return "UNCLASS"

With the expression as:

Reclass(!Mag!,!AP!)

How does this translate into the QGIS Field Pyculator?


texts: QGIS 3 shapefile attribute table sort order

texts:
If I open a shapefile in two empty QGIS projects and view the attribute table the rows are in completely different orders.

How does QGIS decide which order to display the rows in an attribute table and is there a way of making a shapefile display its attribute table in consistent row order? I'm using QGIS 3.4.

Just checked in version 2.14 - row order is consistent so it seems to be a version issue?


texts: Make polygon appear as point when zooming out beyond set range in ArcGIS Online

texts:
Is it possible to have a polygon layer visually appear as a point when zooming out beyond a certain extent that i will determine in a web map in ArcGIS Online?
Or do i have to add a different layer with a point of the center of the polygon and  using the visibility range for each layer have the polygon appear after the point disappears?


texts:
I need to display the heights of the buildings from a given map using QGIS.

I do have OS vector map district data from ordnance survey (https://www.ordnancesurvey.co.uk/) web site that is supposed to contain the heights of the buildings. When I load the data to QGIS, I can see the blocks of the buildings in 2D, but could not find their heights.

I have two attributes, the ID and the feetcode, but the feetcode is the same for all IDs, does it mean that my data does not include heights?

There is actually another file named "TL_SpotHeight.shp", when I open its attribute table I can find an attribute that is called "Height". The problem is the ID and Featcode (it's actually "featcode" non "feetcode", sorry for confusing the words) are all different from the ID and Featcode of TL_Building.shp.


texts: How can I see the layers in an .OSM file?

texts: I have an osm file. I am trying to convert it to geojson with ogr2ogr. I get an error that ogr2ogr can't handle multiple layers. How can I tell what layers are in the file and remove them? I see from How to remove layers of detail from an osm file? that OSM files might not even have layers. So what is the meaning of this error message?


texts: How to verify that all dependency libraries for PostGIS are installed correctly and working correctly?

texts:
I am attempting yet again to build PostGIS from the source to enable the fabled address-standardizer and the reverse-geocoder plugins for PostgreSQL.
I have carefully discovered the source repositories for all the dependencies as listed in the PostGIS manual and I have created local git repositories for each of these.
I will attempt to 'make' them as described on a Linux server machine running ubuntu 18.04. I have already installed PostgreSQL 12.2 from the PostgreSQL.org downloads, along with the server and client tools and additional developer downloads for PostgreSQL.
However, if they build, I would like to also be able to verify they are working as expected. As an experienced developer, I would expect to be able to locate a set of prepared unit tests for this type of requirement.
Does such a suite of 'build tests' exist for the entire group of dependencies listed in the PostGIS manual?
Are there instructions for setting up such tests in the Eclipse CPP implementation?
Are there any workarounds for verifying the operation of dependencies that have no test suites?
Is there a 'standard method' for verifying the operational validity of a CPP project built by one or more make versions that I should be searching for?
I have been using the O'Reilly book service as a source of technical guidance for many years. However, I have been unable to compose a search filter to locate books that possess content that addresses this question.
If any of these questions have 'no' as their answer, or if 'unknown' is the answer, I would be concerned that it will be impossible for me to trust the PostGIS extensions that were compiled according to the instructions in the latest PostGIS 3.02/3 manual.
Is that a reasonable conclusion to make, if I am faced with such an unverifiable situation.


texts: Creating compressed GeoTIFF with QGIS Clipper tool?

texts:
I have a geo-tiff created by Pix4D which is about 375 Mb in file size.

I have clipped it to a smaller area using Qgis 2.8(using Raster ‣ Extraction ‣ Clipper), and the output file is about 900 MB in size.

Why did the clipping process generate such a large file, when the source file was smaller?

How do I run the clipper tool so that it produces a smaller size?


texts:
I would like to detect whenever a user tries zooming but fails due to reaching minZoom/maxZoom of the map. From what I understand, zoom/ zoomstart/zoomend will not be fired as there is no actual animation.

Maybe there's an easy solution. But if not, I could perhaps fire an event manually whenever something happens that could be triggering zoom. I haven't figured out how, though. Probably this is all wrong:

L.Map.ScrollWheelZoom.addInitHook(function() {
      if(map.getZoom() > map.getMaxZoom()) {
        alert("You zoomed too far!");
      }
    });


texts: Change specific tile image source in OpenLayers 5

texts:
I defined a XYZ source and added to map

var source = new ol.source.XYZ({
  url: 'http://blabla.com/{x}/{y}/{z}.png',
  crossOrigin: 'anonymous'
})

During runtime I can access specific tile by

map.getLayers().getArray()[0].getSource().getTile(2,1,-3)

Now I want to change this tile image. I tried setting src_ and image_ properties but no change.


texts: Proper way storing 4D WGS84 routeinformation

texts:
I was wondering if it's possible to store 4D (lat, lon (wgs84), altitude and time) in PostGIS geography type.
I want to store flight paths, where every path consists of multiple points and every point has lat,lon,altitude and time as dimensions. According to the doc I can do this with geometries but I'm wondering what the proper way of doing this would be, since geography seems the right choice over geometry for flight paths that could cross the whole world.


texts: Programatically identifying Joined Field in ArcMap?

texts:
Is it possible to programmatically identify the Join Field that is being used to Table Join two datasets together in ArcMap?   I am currently using ArcGIS 10.0, SP5 and would prefer an ArcPy solution, however I wouldn't be opposed to other solutions, if an ArcPy solution is not available.
One method I tried was looping through all of the fields and look for a matching "baseName", but that is just an "educated guess" where you are hoping that the fieldnames in both databases are the same.
For a graphic representation of what I am after, I basically want to  identify the "Input Join Field" and "Output Join Field" as seen in the "Add Join" dialog, but after-the-fact, of course.

This is a tag-on question to Can a “Join” be detected programmatically?, but in this case I want to extend the functionality to identify the FIELD(s) being used to join the two (or more) datasets together.


texts: "Pagination not supported" error occurs when searching feature layers using the search widget

texts:
I just upgraded to the 3.13 version of the JavaScript API, and I'm working on getting the new search widget (esri.dijit.Search) working in our application.

We have our own ArcGIS Server (10.1 I think) which hosts feature services, and I wanted to be able to search these using the search widget, in addition to addresses and such.

However, every time I try to search one of these feature services, I get a "Pagination not supported" error returned from the server. After a little debugging, I came to find that this was alleviated by removing the "resultRecordCount" parameter of the request that is sent in the query.

I snooped around in the code of the search widget and found that it's running a QueryTask under the hood, and passes a Query with the "num" property set to either the default value of 6 or a value I specify. In the documentation for Query (https://developers.arcgis.com/javascript/jsapi/query-amd.html#num), it's mentioned that the num property is only applicable for services hosted by ESRI.

The workaround right now seems to be creating my own locator service which wraps the feature query functionality, but I'm hoping this issue will get fixed in subsequent API versions. Or is there something I'm missing?

EDIT

Although I can't use any of our feature services in an example, since they're private, suppose we have a reference to an esri.layers.FeatureLayer layer which is hosted on an ArcGIS server version 10.1. What I'm doing is essentially the following:

var sources = searchWidget.get('sources');
sources.push({
    featureLayer: layer,
    searchFields: [layer.displayField],
    displayField: layer.displayField,
    name: layer.name,
    enableSuggestions: false
});
searchWidget.set('sources', sources);
searchWidget.startup();

searchWidget is a reference to an esri.dijit.Search which is placed in the UI. I initialize this widget with the following options:

{
    enableButtonMode: true,
    addLayersFromMap: false,
    enableInfoWindow: false,
    map: <esri.Map reference>
}

Given this setup, my expectation would be that I wouldn't get an error on the web request resulting from performing a search with the layer source selected, but instead I get the "Pagination not supported" error.


texts:
I need to group 400 points in 10 groups with 40 points each. The goal is to get the 40 points geographicly together like this:

Any idea how to do this with QGIS?


texts: Inputting legal description into ArcMap using COGO tool

texts:
I have a legal description that I need to create in ArcMap, and I'm not certain how to input some of the descriptions.

Here are a couple of the descriptions I need to input:

If I only input the part that states:

  "...the chord of which bears North 47° 52' 22" West 88.90 feet to a
  point of reverse curve"

From each section and leave out:

  "...thence 90.89 feet along the arc of a 125.00 radius curve to the
  left"

Then it roughly follows the correct path with a straight line, but it should be slightly curved.

I believe each section is describing a single curved line, but there are 4 parameters, and with the COGO tool, the most I can enter for a curve is three!

So, for the underlined blue and green descriptions in the above image, if I enter these, then it does create straight line segments that roughly follow the path:

Blue:

Green:

But how do I add the

  "...thence 90.89 feet along the arc of a 125.00 radius curve to the
  left"

So that it creates the curved line?


texts: Creating polygon from two not connected linestrings using shapely

texts:
I'd like to know if there is a way to create a polygon in shapely from two not connected linestrings. I've been struggling with this question from some hours now and I can't find a way to deal with it.

For example, I'd like to have the polygon with those linestrings as its exterior coords.
from shapely.geometry import Polygon

p1 = Polygon([1.32, 7.66], [1.4, 5.56], [2.98, 4.5], [5.86, 4], [12.16, 4.88])
p2 = Polygon([1.62, 0], [2, 1], [3.16, 1.68], [5, 2], [8.44, 2.14], [9.92, 0.58], [12, 0], [13.62, 0.8])


texts: Creating a Raster of Road Distance to Point

texts:
I am using QGIS 2.6. This question is in two related parts:

I'm trying to create a raster of travel time from one point to another point travelling via the road system. Each raster pixel will roughly represent the amount of time taken to get to a given point.

Could I do this as a function of the travel-time from a number of different points?  (I suppose I could using the raster calculator once I have the previous rasters.) I.e. I could estimate the total time per week if one travels to location A 5 times a week, location B 3 times a week and location C 1 time a  week.

Additionally, can one incorporate a time for a given distance to the road for points that are off-road (i.e. walking distance)? I.e. an amount per metre of ground to cover the distance to the road (taking the shortest distance possible.


texts: Native algorithms not available in standalone QGIS when 3rd party algorithms are working

texts:
I'm writing some standalone python code (in Spyder) to automate QGIS work. I am successfully using 3rd party algorithms, having performed the following setup:
#initialise qgis
qgishome = 'C:/Program Files/QGIS 3.16/apps/qgis/'
QgsApplication.setPrefixPath(qgishome, True)
app = QgsApplication([], False)
app.initQgis()

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

#import third party processing plugins
sys.path.append(r'C:\Users\weedingb\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins')

from processing_umep.processing_umep_provider import ProcessingUMEPProvider
umep_provider = ProcessingUMEPProvider()
QgsApplication.processingRegistry().addProvider(umep_provider)

However, I can't access 'native' algorithms. For example, running the following (as suggested by the QGIS cookbook)
for alg in QgsApplication.processingRegistry().algorithms():
    if 'buffer' == alg.name():
        print("{}:{} --> {}".format(alg.provider().name(), alg.name(), alg.displayName()))

yields nothing, and inspection of the algorithms list shows none with 'native' as the provider.
Are these not available in standalone python?


texts:
In ArcMap I add 2 tables, join on a field then run an attribute query to return any of the joined tables records that are not null in a particular field - all works OK.

I would like to build this into a model but the model gets stuck when selecting attributes. When I add the same SQL statement "[field] is not null" click verify it says there is an error with the expression but this same expression works when I run manually in ArcMap.

Any ideas why the statement works in ArcMap but not in ModelBuilder Select Tool


texts: QGIS-SAGA Add Raster Values To Points How to apply update

texts:
I have a lattice of points.  I have a surface model.  Both layers share a CRS.  I am trying to populate the point file with the values from the underlying raster using the SAGA tool Add Raster Values to Points in QGIS 2.18 (64bit).  The tool fails with the following message 'module' object has no attribute 'TYPE_RASTER' See log for more details.

I am new to QGIS but not GIS software.

I do not understand why this operation would fail.  I do not understand the error message.

How do I populate my point file with values from the raster?

I see that this is a known bug and there is a fix here but how do I apply the update/fix?


texts:
I would like to store my custom .ui form in the QGIS style saved in PostgreSQL, so I can load it within the layer.

The option is here when you save the style, but I can't see it stored on the database and can' load it.

Has someone already did it? It would help me a lot because I can't install it on each computer.

Maybe a hack by modifying the XML definition in the table?


texts:
I'm working with ArcGIS 10.x and QGIS 1.74 on a Windows XP SP3 System. In ArcGIS I get XML-Data which can be converted to a file-gdb using "IP NAS-Manager".

Included in the XML-files there a owner informations which are stored in a kind of database, labelled filename.data.db, filename.index.db, filename.log.db and so on.

Using ArcGIS 10.x and the extensions "IP ALKIS- Karte" and "IP-ALKIS Buch" it is possible, to show owner information in a separate window. To do this, you must start "IP ALKIS-Buch" before you start ArcGIS 10.x and log on to a Tomcat-Server to verify, that you are a registered user. In ArcGIS 10.x now you select the desired extension, point and click to an item and will get the associated information of owner name, creation-date, length, width etc. of the selected item. All works fine, no problem.

Now I'm trying to get the same information using QGIS, but I don't know the kind of database, where the owner information is stored.

Fortunately, there is another way to get the information. The administration "Kreis" offers "Map Solution" to get the owner information. They also sent me a link to their web-service, where I can log in and search for information. This link can be modified by the item-number I want to search for.

And now comes the problem: How can I realise, that

there is a way to use the link especially to connect to the desired website?
the selected item can be handed over to the link?

The result should be similar to the result of clicking on an item with the "information" button selected and the associated table-information will be displayed.

I hope, I could get you a glimpse of what I'm intending to do, otherwise feel free to ask.


texts: Error running Grass module in Qgis python console

texts:
I am working under Qgis 2.16.0 and i am trying to make a script using the grass module r.water.oulet.
This module has 5 parameters :

1 - Name of input raster
2 - coordinates of outlet point (x,y)
3 - area of the Grass region (xmin, xmax, ymin, ymax)
4 - Celle size of Grass Region
5- Output file

When i use the module in Qgis I don't need to set the area of Grass region. But if d'ont set it in python command :

processing.runalg('grass:r.water.outlet', "D:\Developpement\Debits_cartes_alea\Donnees_test\flow_direction.tif",936400.018,6432568.853,None,0.0,"D:\Developpement\Debits_cartes_alea\Donnees_test\1sortie.tif")

I have this error message : "Error: Wrong parameter value: None"

So I try to set the required coordinates like this :

processing.runalg('grass:r.water.outlet', "D:\Developpement\Debits_cartes_alea\Donnees_test\flow_direction.tif",936400.018,6432568.853,['933963.950608,940746.141038,6431270.66542,6434790.6058'],0.0,"D:\Developpement\Debits_cartes_alea\Donnees_test\1sortie.tif")

Then I have this error message : "Error: Wrong parameter value: ['933963.950608,940746.141038,6431270.66542,6434790.6058']"

Does that come from the syntax or is ther the problem with the values (note that I try with the same values in module and it's work !)


texts: how to identify cross streets for points that are not at intersections?

texts: I have points that fall in the middle of city streets (so not at an intersection). I would like to automate (ideally, with python) a way to identify the cross streets that these points are between. For example, in the screenshot, the point falls on West 231 Street. I would like to find a way to automate a result that returns Palisade Avenue and Independence Avenue. My points file has columns for these cross streets, with the ultimate goal of placing those selected values there.


texts: Calculate length of 2,000 paths - lines in Google Earth?

texts:
I have ~2,000 paths/lines drawn in Google Earth and I would like to know their cumulative length. How can I do this easily?

I have been able to convert into excel the Lat & Long of each point in the string, but some are too long for Excel to accept.
Dylan


texts: What is replacement for setColorShadingAlgorithm in QGIS 2.6?

texts:
I have been struggling with setting the colorshadingalgorithm for a raster layer in QGIS with no success for the past few days. The QGIS cookbook seems to be pretty outdated.

The code that I am using is:

QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()
app = QgsApplication([], True)

plotdir = r"/home/qgis_prac/"
inShp = r"/home/qgis_prac/shps/Peru_Boundary.shp"
uri = r"/home/qgis_prac/tiffs/Peru_DEM.tif"
layer = QgsRasterLayer(uri, "Peru_DEM")
QgsMapLayerRegistry.instance().addMapLayer(layer)
layer.setDrawingStyle('SingleBandPseudoColor') # main thing goes wrong here, I think
layer.ColorShadingAlgorithm(QgsRasterLayer.PseudoColorShader) # seems to do nothng. Also the setColorShadingAlgorithm is deprecated

if hasattr(layer, "setCacheImage"):
  layer.setCacheImage(None)
layer.triggerRepaint()

The code goes on....

My problem is that when the script runs, it gives me an empty map with no entries in the legend. If I use the grayband shader then it works, also for the PalettedColor (but some cells are assigned colors).

What I need is some replacement or work around setting the ColorShadingAlgorithm.


texts: Is it possible to edit layers from a FGDB - File Geodatabse in QGIS?

texts:
I think I've found the answer (you can't edit FGDB's) but there seem to be few reports around that it is possible to edit FGDB files in QGIS.

Can anyone advise if it's possible to setup QGIS (presently 1.8) to edit FGDB's, and how?

It seems like this lack of editing ability could either be a bug in OGR or something to do with having to compile with the FGDB API (I've only ever used OS4GEOW to install QGIS etc).


texts: Displaying null data in choropleth using ArcGIS Desktop?

texts:
I am working in ArcMap and I have created a choropleth map (using a "population" attribute) that is rendering things basically correct except that when the "population" attribute is null, representing no data, it does not render the associated polygon at all.

I would like the null data polygons to render white with a grey outline, as I have the polygons with data rendering with a shade of purple (depending on the "population" value as mentioned above) and a grey outline.

Is there a better way to go about it other than copying the layer and sticking it under the choropleth layer?


texts: Filter and visualize in DataSlider only the dates when images available

texts:
I have a code which is calculating NDVI for Sentinel data and creating DateSlider for this image collection for selected study area. Is there a way to filter and visualize in DataSlider only the dates when Sentinel images available?

https://code.earthengine.google.com/910fdca59ff1a98c9b27cd997d110401


texts: Get a list of property elements using GeoServer

texts:
I'm using GeoServer with openlayers and I want to build a filter check-box dynamically, depending on the elements of this column:

How can I query GeoServer so that I get a list of the different elements in this column?
Or at least all the information of all the columns in XML so that I can parse it?


texts:
I am curious if anyone knows if, and how, it would be possible to configure GeoServer to provide download (GetCoverage) capabilities to multiple adjacent rasters.  A GeoServer mosaic comes close, but provides only the single process image rather than the multiple "unaltered" raw originals that were used to build it.

Ideally, I'd like the user to be able treat the dataset as if it were uniform (e.g. DRGs for a state), but be able to access the original QQ that compose the mosaic.  From what I can tell, GeoServer doesn't support this particular kind of interaction, but perhaps there exists a plugin I haven't yet heard of.


texts: Producing smooth and consistent contour lines from SRTM

texts:
I downloaded and merged together a set of SRTM imagery. I am not sure whether it is my workflow, the actual terrain in question or if the scale of the imagery is simply too large to obtain smooth contours. To further explain, I'm not sure if the area in the image is zoomed out too far thus giving me a ton of small tiny and not connected lines.

For reference, here is the merged DEM raster and my contour results trying different settings:

The tool I am using to run is the "R.contour.step" in QGIS. I set the parameters in various ways but for the most part I left the max contour level at 10,000. Increment at 100 and kept the minimum number of points for a contour line above 300.

If you look at the contour images, there are a lot of small artifacts and lines that do not connect with each other. I am hoping to achieve continual and fluid contours. Running through several iterations has not gotten me any closer to that visualization.

Here is an example of what I'd like to achieve:

Is there a workflow or process I need to go through to achieve the results I am looking for or is this simply a reflection of the terrain in the imagery?


texts: Records in related table do not appear when using the 'Related Tables' button in attribute table

texts:
ArcMap 10 - Geodatabase with related tables & relationship classes.

I've written a script that parses data between two databases. When I go into the attribute table of the point feature class and select a record, I have the option to click the "Related Tables" button.  When I use this button, the resulting table is either empty, or has only 1 record.

The original data has many related records, and they have been parsed correctly into separate feature classes. The records are ALL in the related table, but they're not able to link for some reason. The relationship is through the Date/Time stamp, and I've checked to make sure that the values are the same.

I tried to troubleshoot by adding a join, and everything worked as expected.  Obviously with Joins you can only have a 1-1 relationship as opposed to a 1-Many.


texts: How to find Latitudes and Longitudes of projections of the vertices of a rectangular plane below earth's surface?

texts:
I want to find out the latitudes and longitudes of projections of the vertices of a rectangular plane inside the earth's surface. I know dimensions of rectangle, angles of orientation and latitude and longitude of the projection of the center of the rectangle. Earth's curvature is also to be considered

I added an Image which is a 2D representation of what i want. The plane actually is having 2 angles (its trace makes an angle 'A' with North and the plane makes an angle 'B' with the horizontal). I want to find lat longs of all 4 projections

Note: The dimensions shown are just for depiction. they are not to scale


texts: Editing data in QGIS Server via WFS in AutoCAD Map (2015)

texts:
I have a QGIS Server up and running. I can consult my layers using WMS protocol and edit them using WFS protocol in QGIS or lizmap web client for example, but I am not able to make it work in AutoCAD Map. These are the steps I follow:

In the third screenshot you can see the layer is added to the workspace but the layer is totally empty: it has no geometric data and no alphanumeric data.

I am getting the following errors in the messages window:

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Error al conectar.
Error al abrir la conexión FDO.
Project file error: Error reading the project file

Error al conectar.
Error al abrir la conexión FDO.
Project file error: Error reading the project file

Error al conectar.
Error al abrir la conexión FDO.
Project file error: Error reading the project file

Ha surgido un problema al acceder a la información de esquema de la conexión.
Project file error: Error reading the project file

Error al conectar.
Error al abrir la conexión FDO.
Project file error: Error reading the project file


texts: Is there an equivalent to "keep inner boundaries" in QGIS 3?

texts:
At QGIS 3.4.7, I have some polygons that are overlapped and I use Polygon dissolver (by attribute) from SAGA to dissolve then and I check the option "keep inner boundaries", but this tool takes too much time to run.

I need them that way for the Check Validity to point to the overlapping location.

So, Are there any other way to dissolve the polygons and keep the inner boundaries to gain some time?
The native Dissolver and GDAL´s don´t have that option.

EDIT: Example of inner boundaries in a single feature


texts: Create new polygon layer based on areas of varying numbers of stacked polygons

texts:
I have several species range maps (as polygons) that overlap/stack.  From this, I am trying to identify areas of highest biodiversity.  Ideally, I'd like to create a new polygon layer, where 1 part is the area where 20 ranges overlap, part 2 is the area where 19 ranges overlap, etc.

I am using QGIS and Spatialite as my primary tools.

Thanks for any help that might be offered.


texts: Automatically vectorizing data from WMS using ArcGIS Desktop (ArcScan)?

texts:
I am facing a problem, when I have a historic map in WMS layer and I need to get buildings from this map (approx. 2000) to vector layer fo further analysis.

Is there any workflow how to automatically vectorize these objects from WMS (for example with ArcScan)?

I don´t want to do it all manually, because it is very time consuming.


texts: Python/GDAL: Best way to read in raster for use in A* search algortihm

texts:
I want to create a least cost path program with Python and GDAL since I couldn't find any open source cost path additions for use in Python.

As recommended here the best approach is to use A* search algorithm using slope as the cost between generated nodes.

I want to read in my Slope Raster with GDAL.
What is the best way to read in the Raster in order to further process it in my (not written yet) python A* search algorithm?

So far I read it in as a array but I am not sure if this is the best:

import gdal

# Open data
input_raster = "slope2.tif"
raster = gdal.Open(input_raster)

# read raster as array
bandraster = raster.GetRasterBand(1)
datamask = bandraster.ReadAsArray()


texts: Finding footprint data for US buildings?

texts:
Are there any sources for free and distributable building footprint datasets for major cities around the United States?

Ideally, the data would also have the number of stories of each building; is that available as well?

I'll give you an idea of what I'd like to do: I want to estimate population within an area of a map based on building square-footage and day/night population densities.


texts: Is there a way to calculate the area of a cross section (blue zone on the graph below)in QGIS without creating new layers ? Like a more complex version of the profile tool which would calculate this area from a digitized line?


texts: Totaling polygon area with overlapping buffers in ArcGIS for Desktop?

texts:
I'm using ArcGIS 10.3 to classify landcover within a 30-meter buffer of houses. I categorized this with the Image Classification toolbar, then converted the raster to polygon for calculation purposes. However, some of the houses have overlapping buffers, and after a spatial join & dissolve, one of the overlapping houses/buffer rings will "acquire" any intersecting polygons. This results in overcalculation of the landcover type for one buffer ring and undercalculation for the other (the selected blue polygon in the attached image illustrates this - it attached to the rightmost buffer).

The question is, how do I preserve the integrity of the individual buffers when calculating area for individual landcovers? In other words, how would I go about including the intersected area in each of the buffer rings, while maintaining only the proper total area for each buffer ring (3.14*30^2).


texts: Size of an object in raster image based on DPI resolution and number of pixels

texts:
I'm trying to estimate the size of an object in a raster image, with known DPI resolution.

To evaluate the validity of this method, I do as follows:

1) I'm drawing a 0,02 in × 0,02 in black square in Scribus and export the image to TIFF with 2400 DPI resolution

2) I'm importing the image into Idrisi Selva with GEOTIFF/TIFF import method

3) I'm doing a RECLASS operation to binarize the image, so that square pixels are "1"s, and the background is "0"s

4) I'm running AREA database query operation to estimate the size of the square in "cells". The result is: 1722.

5) Because there are 2400 pixels in one inch, one pixel has 1/2400 inches in length, and thus (1/2400)² area in in²

6) Thus the estimated size of the square is ~0,0012 whereas the real size is 0,0004.

Why the difference? The results are the same when exporting with 1200 DPI, so it doesn't seem to be a matter of accuracy.


texts: AGOL View layer symbology is not updated when updating hosted feature layer

texts:
I have a feature layer called map in AGOL. I want to change the symbology, so I open arcgis Pro and I open this map and overwrite it with new symbols.

I see that the feature layer in AGOL is updated correctly. Only the views of this feature layers aren't updated at all.

Why is this and how can I make sure the views are also updated as well ?

PS: The legend is also not updated of the view layer.

Hosted feature layer (above)

View of hosted feature layer (below)


texts: How to Pass Multiple Dissolve Fields in VBA Code of ibasicgeoprocessor.dissolve Method

texts:
When we use the ibasicgeoprocessor.dissolve method it supports only one dissolve field but if we pass more than one dissolve field then it gives an exception: value does not fall with in the range exception

Please give me any VBA, VB.Net, C#.Net code example of the ibasicgeoprocessor.Dissolve method having multiple dissolve fields.


texts: R: rasterize SpatialPolygonsDataFrame and keep factor field

texts:
I want to rasterize a SpatialPolygonsDataFrame and keep the values of a factor column.

Consider the following code (from Processing vector to raster faster with R):

library('raster')
library('rgdal')
# Load a SpatialPolygonsDataFrame example (Brazil administrative level 2) shapefile
dat <- raster::getData(country = "BRA", level = 2)

# Convert NAMES to factor
dat$NAME_1 <- as.factor(dat$NAME_1)

# Define RasterLayer object
r.raster <- raster()

# Define raster extent
extent(r.raster) <- extent(dat)

# Define pixel size
res(r.raster) <- 0.1

# rasterize
ras <- rasterize(x = dat, y = r.raster, field = "NAME_1")

To access the values from column NAME_1 I tried:

ras_factor <- ratify(ras)
levels(ras_factor)

There are no factor values.

Where do I find the original values from column NAME_1?
How can I rasterize a SpatialPolygonsDataFrame and transfer the factor values to the resulting ras object?


texts: Geometry Generator - make line to polygon perimiter

texts:
I would like to improve the way arrows point to polygons when moving labels. I have some fields with unique numbers associated with each one. When the centroid of a polygon is outside the physical feature, the arrow will be misguiding. How can I improve my expression for the arrow to point to the perimiter of the polygons as illustrated with the blue circles?

CASE
WHEN
("auxiliary_storage_labeling_positionx" IS NOT NULL AND "auxiliary_storage_labeling_positiony" IS NOT NULL) THEN make_line(centroid($geometry),make_point("auxiliary_storage_labeling_positionx","auxiliary_storage_labeling_positiony"))
ELSE
NULL
END

EDIT:

Can anything be done to prevent the arrow appearing a second time if the centroid does not sit on the surface of the polygon as shown below?


texts: Spatial interpolation in PostGIS without outputting Raster

texts:
I have some data that is similar to rainfall, ie. a latitude, longitude, and a Z number. It is in a point layer that is a nice clean point grid. Many Z numbers are missing in this grid.
I want to interpolate to find reasonable Z numbers. TIN, IDW,or other methods. I do not want to go into a Raster first then back out the intersections from the raster onto the point grid.
I can perform an interpolation in QGIS, but it generates a raster which I'm wanting to avoid. The goal is to perform all computation within PostgreSQL for production purposes.
Is this possible?


texts: Visual Studio not creating ".esriAddIn" file

texts: I have an addin that I've been maintaining with Visual Studio 2010 Express for quite some time with no issues. Recently my computer was replaced and Studio 2012 Express was installed. Since that time when I build the addin everything is created except for the .esriAddIn file which is necessary to run the addin with the changes. After looking at the ArcObjects (VB.NET) system requirements I noticed that the only Express version listed was 2010 so I assumed 2012 was the culprit and had our IT department roll me back to 2010. Unfortunately that didn't solve the problem and I have been searching for a solution with no success. Anyone have any ideas why every file is built except the .esriAddIn file?


texts: I do have an overlapping polygon shapefile of different radius. I want to clip a polyline for each fid individually. How is that possible in ArcGIS?


texts: Using arcpy.mapping to make changes to an open mxd

texts:
I want to use arcpy.mapping to make some changes to a mxd.  I am using the Spyder environment for python.  To change the title I can use the following code:

import arcpy

myMap = r'C:\path_to_mxd\myfile.mxd'
mxd = arcpy.mapping.MapDocument(myMap)
mxd.relativePaths = True
mxd.title = 'Title'
arcpy.RefreshActiveView()
mxd.save()

I can't execute mxd.save() on an open mxd file though.  Is there a way to make changes to an open file and view any code I execute on the fly?


texts: Get Lat & Lng from GeoDjango OSMGeoAdmin

texts:
Django OSMGeoAdmin with PointField uses OpenLayers and it's hard to use.

I want to add two new input elements which are lat and lng. Using jQuery, when values of the input elements are changed, I want to change point on the OpenLayers.

How can I do that?

No complete code or instructions are necessary, main idea is enough.


texts: How to "mosaic" rasters and keep attributes

texts:
I have numerous individual 10m Raster files (with both quantitative and qualitative data) for each state and need to combine them into a national dataset.

However, when using a mosaic it seems that the attributes are not kept. I need to keep the attributes that I've already joined so I can conduct further analyses.

Is there an alternative to mosaic, sort of like the merge or append tools for polygons?


texts:
I'm new to this GIS stuff and I've already got the MongoDB extension working, added a database to GeoServer (2.15.1) store and even managed to add a layer and recognizing the document structure (it's a GeoJSON). However, the map is showing as empty. I also added a PostGIS layer just to be sure and that one works like a charm. Here's an image showing the structure GeoServer identifies:

When I Open the layer Preview this happens (code="internalError"Rendering process failed
BasicBSONList can only work with numeric keys not: [type]):

Any Idea why even if I added a layer it's not showing anything? As anyone managed to succesfully connect MongoDB and GeoServer?

EDIT: Here's a the mongoDB JSON example (taken from compass)
type:"Feature"
properties:Object
      id:667881
      Categoría:"Almacenes y bodegas"
      nom_estab:"BODEGA DE LA DIRECCIÓN DE PARQUES Y JARDÍNES DE LA DELEGACIÓN AZCAPOTZ..."
      raz_social:"DELEGACIÓN AZCAPOTZALCO"
      codigo_act:931210
      per_ocu:"11 a 30 personas"
      tipo_vial:"CALLE"
      nom_vial:"PRIMERA PRIVADA RECREO"
      tipo_v_e_1:"CERRADA"
     nom_v_e_1:"SEGUNDA CAMARONES"
      tipo_v_e_2:"CALLE"
      nom_v_e_2:"CAMINO DEL RECREO"
      tipo_v_e_3:"CALLE"
      nom_v_e_3:"TERCERA PRIVADA RECREO"
      numero_ext:656
      letra_ext:null
      edificio:null
      edificio_e:null
      numero_int:null
      letra_int:null
      tipo_asent:"COLONIA"
      nomb_asent:"EL RECREO"
      tipocencom:null
      nom_cencom:null
      num_local:null
      cod_postal:2070
      cve_ent:9
      entidad:"Ciudad de México"
      cve_mun:2
      municipio:"Azcapotzalco"
      cve_loc:1
      localidad:"Azcapotzalco"
      ageb:"472"
      manzana:7
      telefono:null
      correoelec:null
      www:null
      tipounieco:"Fijo"
      latitud:19.47592545
      longitud:-99.18356323
      fecha_alta:"dic-14"
    geometry:Object
type:"Point"
    coordinates:Array
      0:-99.18356323
      1:19.47592545


texts: Reading Esri TIN file directly using QGIS

texts: Can QGIS open directly a TIN file generated with ArcGIS Desktop?


texts: How to call GetLegendGraphics and hide items not in map?

texts:
I need to show a legend from Geoderver layer but wanto to show it only for elements actualy on map. I've found this but seems it not work for me.
Here is my test case:
https://geoservicos.ibge.gov.br/geoserver/ows?service=WMS&REQUEST=GetLegendGraphic&VERSION=1.1.0&FORMAT=image/png&WIDTH=20&HEIGHT=20&LAYER=CGEO:ANMS2010_06_grandesregioes&LEGEND_OPTIONS=hideEmptyRules:true;fontAntiAliasing:true;countMatched:true&bbox=-57.738303574803076,-26.566984791625593,-42.991151777425344,-18.873315675629282&CRS=EPSG:4326

My goal is to hide elements with zero count ( Nordeste, Norte ) in this case.
What I'm doing wrong?
PS: please give me agnostic answers - no OpenLayers / Cesium magic codes... just the URL.


texts: Creating Virtual Layer programmatically in QGIS

texts:
Is it possible to create a Virtual Layer through a Python script? For example, I have a layer "road", and I would like to perform the SQL query
SELECT *
FROM road
WHERE type = 'Expressway'

Will this be possible? Is there any example I can refer to?


texts: Calculate distance from polygon to nearest existing line in geometric network, along network routes

texts:
In ArcMap 10.6, I have a parcels layer and a sewage lines geometric network for a small town. Many of the parcels are still undeveloped. In an effort to accurately assess them for tax purposes, their distance to existing "main" lines affects their value (further from main line = less valuable since connecting parcel to main lines costs a lot).

I need to find distance for each parcel to nearest existing main line along geometric network. In attached image, I need to calculate the distance following the path of the blue line, which follows the route along future mains until it hits existing main. I know that if I did a spatial join of parcels to nearest existing line, it would do so like the red dashed line, which is not what I want.

Is there any way for me to find these distances?


texts: How to get mean values of adjacent cells at point location when extracting values to points?

texts:
I use this code extract raster values to points.

from osgeo import gdal,ogr
import struct

src_filename = '/tmp/test.tif'
shp_filename = '/tmp/test.shp'

src_ds=gdal.Open(src_filename)
gt=src_ds.GetGeoTransform()
rb=src_ds.GetRasterBand(1)

ds=ogr.Open(shp_filename)
lyr=ds.GetLayer()
for feat in lyr:
    geom = feat.GetGeometryRef()
    mx,my=geom.GetX(), geom.GetY()  #coord in map units

    #Convert from map to pixel coordinates.
    #Only works for geotransforms with no rotation.
    px = int((mx - gt[0]) / gt[1]) #x pixel
    py = int((my - gt[3]) / gt[5]) #y pixel

    structval=rb.ReadRaster(px,py,1,1,buf_type=gdal.GDT_UInt16) #Assumes 16 bit int aka 'short'
    intval = struct.unpack('h' , structval) #use the 'short' format code (2 bytes) not int (4 bytes)

    print intval[0] #intval is a tuple, length=1 as we only asked for 1 pixel value

Base on this, I want get mean cells values adjacent at the points. Such as showing in the graph: the value of point is the mean of valid pixels values in the red circle. I assume the radius is n pixels.  Nodata values should be ignored.
Are there any existing function for this processing?


texts: NetworkAnalyst seems to solve along straight lines but not intersections

texts:
I've built a network dataset using a road feature class. All the data I mention is stored in a feature dataset within a file geodatabase.

In particular the connectivity is set to Any Vertex, and not Endpoint. After building the network dataset, I was trying to create an OD matrix using a point feature class. This resulted in no distances other than from a point to itself. This pointed to connectivity problems.

I decided to troubleshoot by manually adding points and solving for a route. Long story short points that are connected by endpoint to endpoint traversal work, but as soon as I add the next point on an edge that intersects, the route solver chokes and states there is no path between the two points.

I'm confused what the problem is since connectivity is set to any vertex. I can verify that it seems the network was built correctly since junctions do exist for both endpoint to endpoint connections as well as intersections.

A similar question is asked here with no clear answer other than a suggestion to make sure that route points be located on the network. I do not believe that is my problem since routes between points are solved in the endpoint case.

The properties dialog on the network dataset says I have:

218035 Junctions
586428 Edges
0 Turns

I'm wondering if the network solver knows it can move from endpoint to endpoint since its just a straight line movement, but without turns for the vertex, there is no way to move onto the intersecting edge.

How can I fix this? I've read about turn feature classes, but am somewhat confused how they apply in this situation. I had thought GlobalTurns would handle the neccesary traversal between any edge provided it was connected.


texts:
Need help with uploading shapefile with srid 102008 to SQL server 2012. I ran the ogr2ogr script in the ogr shell and it uploaded the shapefile within seconds but when I try using python subprocess to call ogr, I get an empty spatial_ref_sys and geometry column table added to SQL server 2012 database. I should add that I have read write access to the database .

Here's my script-

import os
import sys
import subprocess

# Set OGR2OGR environment
gdal = r'C:\Program Files (x86)\GDAL'
gdalData = r'C:\Program Files (x86)\GDAL\gdal-data'
os.environ['PATH'] = gdal + ';' + os.environ['PATH']
os.environ['GDAL_DATA'] = gdalData

indriver = "MSSQL:server=<server>;database=<database>;trusted_connection=yes"
ogr2ogr = r'C:\Program Files (x86)\GDAL\ogr2ogr.exe'
cmd = '-f MSSQLSpatial'
proj_cmd = '-a_srs'
proj = "EPSG:102008"
reName_cmd = '-nln'
reName = 'NewName'
shp = path = r'H:\NA\20151226_LNAD83.shp'
translateFile_tx =' '.join([ogr2ogr, cmd, indriver, shp])
subprocess.call(translateFile_tx)
print 'Done'


texts: Wrapping legend text in QGIS MapComposer at character count

texts:
I want to wrap legend text after a certain amount of characters. Is there an automatic way to do that?
I don't want to enter a character (e.g. *) for each line as described in this post, that is too much work...

Edit: I tried to use an expression, but the button is greyed out...


texts: change ggmap to sp.: Create circle where radius represents distance (km)

texts:
I am trying to made multiple concentric circles of multiple radii in kms scaled to a map, similar this post was using the post here (ggmap: create circles).

map.sp<-readShapeSpatial("C:file.shp")

plot(map.sp,xlim=c(143.4, 144), ylim=c(-13.4, -13))#"SpatialPolygonsDataFrame"

d8 <- data.frame(lat = c(-13.177047), lon = c(143.949718))
coordinates(d8) <- ~ lon + lat
projection(d8) <- "+init=epsg:4326"
d_mrc8 <- spTransform(d8, CRS = CRS(projection(map.sp)))

library(scales) # for `alpha()` function
library(rgeos) #for gBuffer function
d_mrc_bff8 <- gBuffer(d_mrc8, width = 100) # Buffer creation

plot(d_mrc_bff8, col = alpha("blue", .35), add = TRUE)
points(d_mrc8, cex = 2, color="red",pch = 20)

All commands finish without any error messages.

My map is produced within the correct x and y limits and the final point is plotted (although in black not red).

However, no buffer is created.

Further information below about the specific "SpatialPolygonsDataFrame" file I am using.

coord. ref: +proj=longlat +ellps=WGS84

Extent: 138.0012, 153.5516, -29.17855, -10.34855  (xmin, xmax, ymin, ymax)


texts: ST_MakePoint or ST_PointFromText to generate points

texts:
I'm new to postGIS and am having trouble finding an answer to this question. I have loaded UK postcode data (OSGB36, SRID 27700 with eastings and northings fields) into a postGIS database. I now want to generate a geom column to store the point information.

I can find some clear instructions on how to do this using ST_MakePoint from here: http://twiav-tt.blogspot.co.uk/2012/07/postgis-using-latitude-and-longitude-to.html

This gives the nice straightforward syntax example of:

UPDATE netherlands.airports SET geom = ST_SetSRID(ST_MakePoint(longitude,latitude),4326);

I've read that although ST_MakePoint is fast, it's not as standards compliant/transferable as ST_PointFromText.

But no matter where I look for a similar example for ST_PointFromText, I can't find one. For example, on Boston GIS, the following example is given:

UPDATE points_of_interest SET thepoint_lonlat = PointFromText('POINT(' || longitude || ' ' || latitude || ')',4326)

The postGIS page gives a near identical example:

SELECT ST_PointFromText('POINT(-71.064544 42.28787)', 4326);

These previous two example just seem to show you how to generate a point for some text that you've hand-typed in, they don't show you the syntax for how to derive it from other fields in your table.

So I know I'm being a noob here, but I have one column called 'easting' and one column called 'northing'. What is the syntax to create my geom column using ST_PointFromText?


texts: "Dataset xxx does not exist or is not supported" using Table to Domain in Python

texts:
I have a directory of hundreds of *.csv files, containing code/description pairs for geodatabase domains. I'm trying to automate the creation of these domains using Python.

If I manually run the Table to Domain tool using a *.csv file in
ArcToolbox, the domain is correctly created
If I create a Model using the Table to Domain tool, hard-coding in the input file (and
all other parameters) then export this to a Python script, the
domain is correctly created when I run the Python script
If I try to iterate through all *.csv files in the directory, running the Table
to Domain tool for each file, I get the ""Dataset xxx does not exist
or is not supported" error message:

  Failed to execute. Parameters are not valid.ERROR 000732: Input Table:
  Dataset C:\temp\01-01.csv does not exist or is not supported ERROR
  000308: Invalid field type Failed to execute (TableToDomain).

I suspect the problem is due to the path not being correct for the *.csv file, but I can't work out how to get it working. os.path.exists shows that the file does exist, so it may be an internal ArcGIS thing.

Note that the *.csv file does work correctly via ArcToolbox and the hard-coded Python script, so it's unlikely that the field type really is invalid.


texts: Build tile server manually with PostgreSQL db on a different server

texts:
I have set up my map tile server and have used the following components:

Mod_tile
osm2pgsql
postgres
Mapnik
open street map carto

In my configuration, postgres is on a different server. when I access the /osm_tiles/0/0/0.png, I get a blank square. I think the data is not being fetched from postgres database. I have mentioned the configuration in the osm2pgsql/project.yaml file.

How do I tell my server that the PostgreSQL db is on another server? Do I need to mention the db details in some other place as well?
Which of tool fetches the tiles from the PostgreSQL db to be displayed?

P.s I've used this link as the setup guide.


texts: Why does a layer imported from GPS have lat/lng rather than project coordinates?

texts:
I am importing tracks from my GPS using the GPS Tool in QGIS. I check the layer CRS and it shows WGS 84 but when I try and view the layer in my project which has a CRS of EPSG:2193 the coordinates are still lat and log (i.e. native WGS84).

I have the same data on another computer and that works fine.  I have compared all the setting that I can think of (Project Properties and Setting->Options) between the two system and all the ones I looked at are the same so far as I can tell.

In particular I have made sure that "settings->options->CRS->'on the fly projection' is checked.

I am pretty sure that I have had this issue before but can't remember how I solved it.

What setting do I need to get QGIS to automatically convert the wgs84 layer to project projection?  Alternatively what setting might I have inadvertently turned on that might have broken this?

later:
I worked through the various things that people suggested (thanks very for the suggestions BTW) I and still could not figure out what the problem was so I then started from scratch again and now everything works as expected.  Sigh...

I have been in IT long enough to know that when things just don't work as they should cycling power (or the software equivalent there of) should be tried before shouting "HELP!"

not sure of the etiquette here -- should I mark this as answered?


texts: Creating time-aware layer (ArcGIS API for JavaScript + ruby rails) using data from server?

texts:
I have a data set of points.
Each point has latitude, longitude, date and description
sent from the server

I already have a graphic layer which is used for represening the points in a static way with some popups
like this

          ...
          var graphic = new esri.Graphic(point, symbol);

          graphic.setInfoTemplate(new esri.InfoTemplate("", "Number: " + data.number));

          gLayer.add(graphic);
          ...

But Now I need to add a time-aware layer (a feature layer I believe) that allows me to show the data according to datetime information using a time slider, like here [oil Filed][1].
The only difference, I need to use the data sent from the server as json
and I need to add popups to each point

I believe I need to add a featureLayer (MODE_SNAPSHOT) with a featureCollection as params that contains the data sent from server
but I coulnd't figgure out how to do this (I'm a newbie in ArcGis, several hours only).

To be exact I'm not sure how to create that feature collection with latitude/longitude for each point and the most important thing: how to connect it with corresponding time information

I.e. I have sets of

{lat: xxx,
long: yyy,
datetime: yyyy-mm-dd,
description: sometext
}

sent from the server

Now I need to figure out how to pass the data to a featureLayer so it can be shown on a map with a time slider  [1]: http://help.arcgis.com/en/webapi/javascript/arcgis/help/jssamples_start.htm#jssamples/time_slider_layerDef.html "oil Field"


texts: Create Random Raster (Data Management) Generating "No Data" values

texts:
I would like to generate a random raster with a cell size of 10 to match an elevation raster's cell size. However attempting to do so with a defined extent creates this incomplete raster with a section of no data:

If I up the cell size to 20 I get a random raster as expected without sections of "No Data":

What is causing this behaviour and how can I prevent it while maintining the matching cell size/raster resolution?

I am using ArcGIS Pro 1.4

The "Create Random Raster (Data Management)" tool was used here, using a Uniform distribution with a range of 0 to 10. The extent of the output was set to the extent of the elevation raster using the drop down menu - this defaulted the cell size box to 72, which I manually changed to 10 to match the actual cell size of the elevation raster (taken from the properties dialogue), perhaps this is a clue:


texts:
I would like to create a risk map using RFE Data and Tiffs but now I realised that I actually have no idea how.

I have dowloaded .Tiff(Region), .hdr and .bil(Rainfall/RFE) files from the USGS page, after that I georeferenced the .tiff and the .bil data to see in which part of the Region there is the most Rainfall( It's about malaria which has the highest risk in regions with the highest RFE.)

What I now want to do is to create a map that is clearly showing the regions where the risk is at its highest and a way to visualize it as good as possible.