texts: Convert a SQL Server Geometry into PostGIS spatial data with Python and shapely

texts:
I'm trying to migrate spatial data from a SQL Server database to a PostGIS database with Python. I need this to be automated since its not just one table, but a lot tables in a lot of databases.

I am struggling with the geometry types.

From the SQL Server I get this:

b'\xe8d\x00\x00\x01\x04\x05\x00\x00\x00\x05\xa4\x1e\xbe\xee\x15\x12Aa/\xb9y)\x97UA\xdc\xd2\xe2\xd6\xeb\x15\x12A\xc73\xd7\xf7)\x97UA \xa8\xb07\xef\x15\x12Ac\xa9\xbf\x0b*\x97UAIy\xec\x1e\xf2\x15\x12A\xfd\xa4\xa1\x8d)\x97UA\x05\xa4\x1e\xbe\xee\x15\x12Aa/\xb9y)\x97UA\x01\x00\x00\x00\x02\x00\x00\x00\x00\x01\x00\x00\x00\xff\xff\xff\xff\x00\x00\x00\x00\x03'

Which according to my research should be a WKB (Well Known Binary)

Inserting this directly into my PostGIS database gives me this error :

ERROR:  Invalid endian flag value encountered.

Which, again according to some research, I need to parse either into WKT or a string before being able to insert it.

Right now my idea is, to use shapely to parse it into a simple string representation and then insert it, but i cannot even parse it, code looks like this :

def WKBtoStr(self, WKB):
    print('befpre')
    print(WKB) # gives the above byte sequence
    geomString = wkb.loads(WKB) # crashes see error below
    print('after')
    print(geomString)

This gives me the following error :

ParseException: Unknown WKB type 100
Could not create geometry because of errors while reading input.

I have also tried :
Decoding the byte sequence, which just yields wierd gibberish.

Any suggestions?


texts:
I´m trying to merge several partially overlapping rasters using Gdal Merge in QGIS 3.14.15-Pi. The merge process runs without throwing any errors but the output is incorrect. The output has no defined CRS, extent and the width and height are incorrect. Basically the output fails (cannot even view it). Why is this happening?
I´ve attached a few images of the same.


texts: Merging 2 non overlapping points datasets

texts:
I have 2 dataset of points, each one measures a different variable, i need to "merge"(not sure if it is the correct term) these two into one dataset.
My first instinct is to take a purple point find k closets green points and assign the average as the green dataset value for the purple point.
Example image of the points.

Is there a recommend method for this kind of problem?
Edit:
I am using R, but if there is not package that solves the problem i will be fine with a theoretical reference.
Edit 2:
Added code to generate an example, i would like a to know which methods are available and which ones are recommended for a case like this.
library(sf)

N_points <-  100
X_range <- c(0, 20)
Y_range <- c(0, 20)

set.seed(123)
df_1 <- data.frame(
x = truncnorm::rtruncnorm(n = N_points,  a=X_range[1], b=X_range[2], mean=10, sd=5),
y = truncnorm::rtruncnorm(n = N_points,  a=Y_range[1], b=Y_range[2], mean=10, sd=5),
value = rnorm(n = N_points, mean=7, sd=2)
)  %>%  st_as_sf(coords=c("x", "y"))

set.seed(753)
df_2 <- data.frame(
x = truncnorm::rtruncnorm(n = N_points,  a=X_range[1], b=X_range[2], mean=10, sd=5),
y = truncnorm::rtruncnorm(n = N_points,  a=Y_range[1], b=Y_range[2], mean=10, sd=5),
value = rnorm(n = N_points, mean=10, sd=2)
)  %>%  st_as_sf(coords=c("x", "y"))

plot(st_geometry(df_1), col="purple", axes = TRUE)
plot(st_geometry(df_2), col="green", add=T)

Edit 3, my description of the problem was not clear enough:
Example:
Each purple point represent a measure of air pollution.
Each green point represent a measure of radiation.
I want as result a datasets that has both a measure of radiation and a measure of air pollution for each point.


texts: Mapquest routing using multiple stops(via) in OL3

texts:
I am trying to find the route between two stops in Mapquest by passing the lat long details in the URL,

http://open.mapquestapi.com/directions/v2/route?key=YOUR_KEY_HERE&outFormat=json&routeType=fastest&timeType=1&enhancedNarrative=false&shapeFormat=raw&generalize=0&locale=en_US&unit=m&from=38.8941,-77.07556&to=38.84458,-77.07823

But I want to find the route that passes via specified location,How to achieve this?


texts:
I have 300 .tiff files, that i need to create raster stack of them. Can you please help me with R code for that?

I can do in Arcgis, Envi, R. I have tried:

setwd("F:\\MODIS\\Modis EVI\\HDF8 EVI")

list.files(path="F:\\MODIS\\Modis EVI\\HDF8 EVI", pattern = ".tif")

tifs <- list.files(path="F:\\MODIS\\Modis EVI\\HDF8 EVI", pattern = ".tif")
stackSave(tifs,"hdf8_EVI.TIF")

and

setwd("F:\\MODIS\\Modis EVI\\HDF8 EVI")
raster_stack=list.files(getwd(), pattern="tif$", full.names=FALSE)
for(i in raster_stack) { assign(unlist(strsplit(i, "[.]"))[1], raster(i)) }

Both didn't work!


texts: Why is only one chart showing in one polygon in ArcGIS (10.3) and not all six polygons?

texts: I have a small dataset for six polygons (regions of the UK) with nine fields of numeric data. I'd like to show the data as a pie or bar chart in each polygon, but having set this up under symbology ArcGIS is only displaying one chart in one of the polygons rather than all six. The polygon it's working for is the first row in the Attribute Table. I can't work out why the other five aren't displaying.


texts:
I have a topographic map with Buildings, Roads, Hydro and land use as the prominent layers in polygon format. Also i have newly collected Building polygons of the same area in shape file format. Now i need to have an overlay analysis of Buildings in these two resources. I can rasterize and geocode the map to orient with Digital data. Also these two can be overlaid and compared manually. But i need to do this analysis automatically with possibly using a shorter method. One possibility is to vectorize Buildings in the map and perform vector overlay.
However I am looking at a more efficient method possibly filtering Buildings in raster and comparing it with Buildings converted to rater in shape file. Any idea on this is appreciated.


texts: Is the bilinear resampling technique appropriate for the projection of ASTER data?

texts:
I need to project some ASTER data from a geographic coordinate system to a projected one. I have read online that the cubic resampling technique is the most appropriate for this purpose. However, using this technique results in the output having negative values, which I cannot have. ArcGIS help instructs the user not to use the nearest, or majority techniques for the projection of continuous data. This leaves only the bilinear resampling technique, which does not have the negative value problem associated with the cubic technique. So, is the Bilinear resampling technique appropriate in this case?

EDIT: I have said that I am using this ASTER data to map a lithic formation which is composed of chert. I am also attempting to use it to discern between dolomitic formations which alternate betweeen chert rich and chert poor.


texts: PyQGIS 3: Create a layer with a primary key

texts: I have  written a plugin for QGIS that is importing some data. In this plugin, I can create a memory layer for the data like I have described here https://gis.stackexchange.com/a/349354/10123, but how can I create a primary key / feature id column on this table? In the QgsVectorDataProvider, I see I can use pkAttributeIndexes to list the fields making up the primary key, but how can I set it?


texts: Is there any Google Maps-like zoomable world map, but with an equal-area projection?

texts:
Whenever I look at Google Maps at a global or polar area, the projection hurts.  Is there any web service with a zoomable world map, that uses an equal-area projection — preferably pseudorectangular such as the Mollweide projection?  Even National Geographics online map uses Mercator.

Another alternative, Google Earth, does not particularly distorts the Earth, but does not give a global view either.

See also: Why has Web Mercator (auxiliary sphere) become the web map standard?

Edit 2018-08-24: Google Maps no longer uses Web Mercator, but has moved to what appears to be an orthographic projection.


texts: lidR: Error in .bboxCoords(coords) : nrow(coords) > 0 is not TRUE

texts:
I'm using lidR to create some lidar metric rasters that match the layout of Landsat images. In the lidR documentation it says this can be done simply by using a RasterLayer instead of a resolution in the grid functions.

las <- readLAS(here("data/lumber1.las"))
evi <- raster(here("data/ls8_evi.tif"))
projection(las) <- projection(evi)
point_density <- grid_density(las, evi)

The code is simple but I get

Error in .bboxCoords(coords) : nrow(coords) > 0 is not TRUE

when I run the grid_density() function and I'm not sure how to get around this. If I use a number as the resolution instead of a raster all the grid functions work.

> print(evi)
class      : RasterLayer
dimensions : 10000, 10000, 1e+08  (nrow, ncol, ncell)
resolution : 30, 30  (x, y)
extent     : 1484415, 1784415, 1364805, 1664805  (xmin, xmax, ymin, ymax)
crs        : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : C:/Users/.../R/lidR_test/data/ls8_evi.tif
names      : ls8_evi
values     : -2.496704, 2.983905  (min, max)

> print(las)
class        : LAS (LASF v1.4)
point format : 6
memory       : 994.8 Mb
extent       :2e+06, 2005000, 230000, 235000 (xmin, xmax, ymin, ymax)
coord. ref.  : NA
area         : 25 kunits²
points       : 11.34 million points
density      : 0.45 points/units²
names        : X Y Z gpstime Intensity ReturnNumber NumberOfReturns ScanDirectionFlag EdgeOfFlightline Classification ScannerChannel Synthetic_flag Keypoint_flag Withheld_flag Overlap_flag ScanAngle UserData PointSourceID


texts: Using different types of barrier lines in ArcGIS Network Analyst?

texts:
I'm trying to load line barriers into network dataset. The shapefile I prepared to use as barrier has a column called "TTtype" with two different types of "Restriction" and "Scaled Cost" and also a column called "TTFactor" with the factor for scaled cost type. When I load locations, the results are only the default values. I expected to see the number of restriction and scaled cost under the Line Barriers in network analyst panel.

What's the wrong with my approach?


texts: Creating SQL Parametric View based Layer in QGIS?

texts:
How to create SQL Parametric View based Layer in QGIS (like in GeoServer documentation)?
I have PostGIS spatial data. There I have my own function, that does something. I execute it from GeoServer:
select * from my_scheme.get_geom_all('%T_KEY%','%KPI_NAME%','%FLTR%')

This function gets params via a web app (&viewparams=T_KEY:2018-12-01;KPI_NAME:any_KPI;FLTR:any_value_1>=FLTR:any_value_2<FLTR:any_value_3 )
On GeoServer it work fine but I have not been able to duplicate it on QGIS.


texts: Let's say I have a print composer with two maps. On one I define a grid. Can I copy/paste this grid over to the other map in any way?


texts: Does anyone know where I could find a class diagram for the different modules of geoexplorer (geoext and gxp)?


texts: Read SRTM HDG elevation files with C++ / Strange negative numbers (not Null)

texts:
I'm trying to write a little C++/OpenGL utility with which I could manage SRTM HDG files. For now, I'm just trying to read the files correctly. So I downloaded a sample tile from the main server, which is associated with this image.

To generate an image like that is exactly what I want to do initially, by code, no lib (if it doesn't became too hard). The problem is that the best I got so far are images like this (compare!).

I have tried a couple of proceedings... Basically, I'm transforming altitude in pixel color (from black/low to white/high), and putting the points per pixel with openGL as if creating an image. I've done it directly already, and now I'm averaging with the points around the point of interest. The results are very alike in both cases.

Before that, when I was visualizing the raw data (the altitude int's), I noticed a lot of strange negative values (like, -100 aside to 512 and so on), and it seems they correspond to some of the dark areas in my images (where they should not be so, if compared to the original image from NASA).

Here is my code for the draw() function (ok, it's just to get something out). Im usign netBeans IDE, on linux/Ubuntu:

void draw()
    {
       glClear( GL_COLOR_BUFFER_BIT );
       int rows = 0;
       int cols = 0;
   byte = new char[2];
   file.seekg (0, ios::beg);
   file.read(byte, 2);
   while(!file.eof())
   {
      alt = byte[1] | ( byte[0] << 8 );
      if (alt < minimo) minimo = alt;
      if ( alt > maximo ) maximo = alt;
      h[cols][rows] = alt;
      cols++;
      if ( cols > 3600 )
      {
         rows++;
         cols = 0;
      }
   file.read(byte, 2);
   }
   for (rows = 1; rows < 3599; rows++)
   {
      for (cols = 1; cols < 3599; cols++)
      {
         altitude[cols][rows] = 0.1111111111*( h[cols][rows] + h[cols-1][rows-1] + h[cols][rows-1] + h[cols+1][rows-1] + h[cols-1][rows] + h[cols+1][rows] + h[cols-1][rows+1] + h[cols][rows+1] + h[cols+1][rows+1]);
      }
   }
   glPointSize(1.0f);
   glBegin(GL_POINTS);
   for (rows = 1; rows < 3600; rows++)
   {
      for (cols = 1; cols < 3600; cols++)
      {
      color = 0.0005555555555556f * (float)altitude[cols][rows];
      glColor3f(color, color, color);
      glVertex3f(cols, rows, 0);
      }
   }
   glEnd();
   cout << "lemos tuo - cols = " << cols << ", rows = " << rows << " - minimo: " << minimo << ", maximo: " << maximo << endl;
   glutSwapBuffers();
}

You may notice, also, the "min" (minimo) and "max" (maximo) values. When I compare with the min and max I get from GRASS, for instance, the max's are the same (1644 m), but my min = -128 (negative), and GRASS's min = 120 (positive)! And GRASS renders it just like SRTM's original.

Can anyway help me? What am I doing wrong? Am I reading wrong or rendering wrong, or both? I've gone through a some official docs and other questions, and all I got is that I must swap the bytes. Some people talk about negative numbers (that should be small, and correctly positioned - there is a river in this image). But I can't manage to get help about some - 100 at the side of some +500's that build strange holes in the image - but GRASS get's it right. So... My mistake. But where?

PS: first time I do things with these kind of dada.


texts: Parsing error: in Conditional raster calculator?

texts:
I was trying to calculate Soil types based on the FAO's guidline using the Con tool in raster calculation with ArcGIS 10.3. from three layers (silt, clay & sand). But, I got an error saying:

  Python syntax error:Parsing error SyntaxError: invalid syntax line1

I have put the code below. What is wrong with it?

OutRas = Con(("siltwtdavg" + 1.5*"claywtdavg")<15,1, Con(( "siltwtdavg"+ 1.5*"claywtdavg" >= 15) & (silt + 2*"claywtdavg" < 30),2, Con(("claywtdavg" >= 7 & "claywtdavg"< 20) & ("sand_av" > 52) & (("siltwtdavg" + 2*"claywtdavg") >= 30) | ("claywtdavg" < 7 & "siltwtdavg" < 50 & ("siltwtdavg"+2*"claywtdavg")>=30), 3, Con(("claywtdavg" >= 7 & "claywtdavg" < 27) & ("siltwtdavg" >= 28 & "siltwtdavg" < 50) & ("sand_av" <= 52),4, Con(("siltwtdavg" >= 50 & ("claywtdavg" >= 12 && "claywtdavg" < 27)) | (("siltwtdavg" >= 50 & "siltwtdavg" < 80) & "claywtdavg" < 12),5, Con("siltwtdavg" >= 80 & "claywtdavg" < 12, 6, Con(("claywtdavg" >= 20 && "claywtdavg" < 35) && ("siltwtdavg" < 28) && ("sand_av" > 45), 7, Con(("claywtdavg" >= 27 && "claywtdavg" < 40) && ("sand_av" > 20 && "sand_av" <= 45), 8, Con(("claywtdavg" >= 27 && "claywtdavg" < 40) && ("sand_av"  <= 20), 9, Con("claywtdavg" >= 35 && "sand_av" > 45), 10, Con("claywtdavg" >= 40 && "siltwtdavg" >= 40), 11, 12 )))))))))))


texts: Georeferencer Plugin: Installed but does not appear in raster menu (QGIS 10.3, Windows 10)

texts: I see that others have had problems finding a rouge georeferencer window, or using the plugin. But when my student tries to use the confirmed installed georeferencer plugin in QGIS 3.10 on Windows 10 no plugin appears in the raster menu. Any thoughts?


texts: Local variance image in python using gdal and a running window approach

texts:
I want a local variance image with a 3x3 of a geospatial raster image using python. My approach so far was to read in the raster band as an array, then using matrix notation to run a moving window and write the array into a new raster image. This approach worked well for a high pass filter as described in this tutorial: http://www.gis.usu.edu/~chrisg/python/2009/lectures/ospy_slides6.pdf

Then I tried to calculate the variance with several approaches, the last one using numpy (as np), but I just get a gray image with the same value everywhere. I am open to any kind of solution. If it gives me the average local variance in the end, that would be even better Thanks, Mace

  rows = srcDS.RasterYSize
  #read in as array
  data = srcBand.ReadAsArray(0,0, cols, rows).astype(np.int)

  #calculate the variance for a 3x3 window
  outVariance = np.zeros((rows, cols), np.float)
  outVariance[1:rows-1,1:cols-1] = np.var([(data[0:rows-2,0:cols-2]),
    (data[0:rows-2,1:cols-1]),
    (data[0:rows-2,2:cols]  ),
    (data[1:rows-1,0:cols-2]),
    (data[1:rows-1,1:cols-1]),
    (data[1:rows-1,2:cols]  ),
    (data[2:rows,0:cols-2]  ),
    (data[2:rows,1:cols-1]  ),
    (data[2:rows,2:cols]    )])
  #output
  outDS = driver.Create(outFN, cols, rows, 1, GDT_Float32)
  outDS.SetGeoTransform(srcDS.GetGeoTransform())
  outDS.SetProjection(srcDS.GetProjection())
  outBand = outDS.GetRasterBand(1)
  outBand.WriteArray(outVariance,0,0)
  ...


texts: Python Script examples for geoprocessing shapefiles without using arcpy

texts:
I would like to use a Python script that is not based on arcpy to do things like  query a shapefile by attributes, create new layer from selection, and calculate areas of a polygon and convert polygons to points.

Anyone have any code examples of using other Python modules or libraries? I am able to do this easily using arcpy but i wanted to explore other options.


texts: Add and Updates not working for attributes in ArcGIS Runtime for .Net

texts:
Environment:

ArcGIS Runtime for .NET 10.2.7
ArcGIS Server 10.3.1
Map Service with Feature Access enabled (Create, Read, Update)
Postgres database underneath

ServiceFeatureTable reportLayer = await ServiceFeatureTable.OpenAsync(new Uri("http://localhost/arcgis/rest/services/MyMapService/FeatureServer/4"));
var testFeature = reportLayer.CreateNew(); //new GeodatabaseFeature(reportLayer.Schema);
testFeature.Attributes["a_commodity_group"] = "fake";

var addedReportOIDs = await reportLayer.AddAsync(testFeature);
var editResult = await reportLayer.ApplyEditsAsync();

I've done the CreateNew() as well as creating a new GeodatabaseFeature as described in the documents here:

https://developers.arcgis.com/net/desktop/guide/edit-features.htm#ESRI_SECTION1_2A4C0222216A448C9CBAFCD0F75A21F6

I've tried Adding and Updating features. The "a_commodity_group" attribute is exactly how I see it spelled in the rest endpoint and is a Text/varchar field in the database underneath

I'm seeing none of my attributes come through. (example, the "fake" value that I've put in here)

ArcGIS Server Log (at "Info" level)

  FINE Mar 28, 2016, 8:59:46 PM GraphicFeatureServer::HandleRESTRequest

  inputJSON = {"rollbackOnFailure":true,"adds":[{"attributes":{"objectid":-10}}]}

Inspecting client traffic from Telerik's Fiddler:

As you can see, the request is totally wrong in that it doesn't have any other attributes going in.

Manually adding features through "ApplyEdits" in the REST endpoint worked.

Info from ArcGIS Server logs (the attributes come in correctly here):

  GraphicFeatureServer::HandleRESTRequest ## inputJSON =
  {"adds":[{"attributes":{"objectid":405,"a_commodity_group":"fake"}}],"rollbackOnFailure":true}

I don't see any known issues with this in the ArcGIS Runtime.
I'm considering writing my own REST api as the workaround right now since I can't get this to work.

It's quite possible that I'm making a totally wrong assumption in how I'm using the API, but I haven't seen many examples for adding features to a connected service. I'll double check the full .net runtime sdk pdf document and see if I'm doing anything wrong, too.

Update:
I also noticed that DeleteAsync doesn't work unless I layer.QueryAsync the features first. Probably has something to do with some local cache in the ServiceFeatureTable layer object. Before I was querying objectids using the ExecuteObjectIdsQueryAsync() method on a QueryTask object.


texts: Auto populating attribute table with lat and lon in QGIS

texts:
is there a simple way (without writing script or Python skills) to auto calculate and auto populate lat and lon attributes in QGIS?

When I manually add new points or move existing points in my layer, I'd like the lat and lon attributes to be captured and the attribute table populated/updated as I complete the task. Right now, I create or modify the points and then update the attribute table with $y and $x in the Field Calculator.


texts: How to change property of features in a FeatureCollection where there is overlapping?

texts:
I have two feature collections and I want to change the features' property in one collection based on the other one where there is overlapping, but it seems that I just merge and keep both features where they overlap.
Map.addLayer(Grid,{color:'red'},'Grid')
var Fcorrectioend=function(feature){
                    var day=ee.Number.parse((ee.String(feature.get('event_date'))).slice(8))
                    var Month=ee.Number.parse((ee.String(feature.get('event_date'))).slice(5,7))
                    var Year=ee.Number.parse((ee.String(feature.get('event_date'))).slice(0,4))
                    var date=ee.Date.fromYMD(Year,Month,day)
                            return feature.set('system:time_start',date)
                              };
var test=table.map(Fcorrectioend)

var day1=test.filterDate('2020-01-01','2020-01-02')
Map.addLayer(day1,{},'day1')
print(day1)

var test=(day1.merge(Grid)).select('event_count')
Map.addLayer(test,{color:'green'},'testMerge')

Here is the link to the code:
https://code.earthengine.google.com/7378c612af5339f972f6fc0b0e9ad85c


texts: Authagraph world map with deformed south pole instead

texts:
Read this first:
The north pole is deformed on AuthaGraph world map

I am interested in an authagraph map, where the tip of the tetrahedron is orientated to the south pole instead, so the deformation would be condensed there.

Did someone generated such a map?

Or if not how could i generate one?


texts: Add and later access Leaflet.draw polyline index and a new property value

texts:
I am using the Leaflet draw plugin located at

https://github.com/Leaflet/Leaflet.draw

and what I want to do is when a person draws a new polyline to print in the console the index number of that polyline and also add a value such as color:'red' to its properties.

Then when I click that line to edit it later also print to the console the index number and the value stored in the color property.


texts: Searching for value in layer without opening MXD file?

texts:
As part of the validation for parameters of a Python toolbox, I need to check if the text in parameter #1 is exists in a layer in a MXD file selected as parameter #2. Currently, I use this code in updateMessages:

mxd = arcpy.mapping.MapDocument(parameters[2].valueAsText)
layer = arcpy.mapping.ListLayers(mxd, "LayerName")[0]
condition = u"Field=" + parameters[1].text
with arcpy.da.SearchCursor(layer.dataSource, ["Field"], condition) as cursor:
    #Do some stuff...

The problem is that this takes very long time since the rather large MXD file has to be opened. It is not acceptable to have a two minute waiting time with ArcMap freezing every time the parameter validation is performed.

So my question is this: Is there any way to find the source of the layer, and search for a value in it, without having to open the MXD?


texts: Populating attribute table string field with values from long field using vbScript in ArcGIS Field Calculator?

texts:
I want to calculate the field "area_type" with text depending on numbers from the "Value" field. I have tried with the following VB-script in "Field Calculator" but am getting this error message:

"The calculated value is invalid for the row with ObjectID = 1. For example, the calculated value may be too large for the field or you may be trying to add a string to a number field. This row will not be updated".

My "area_type" field type is string and "Value" is type long.

Dim area_type

Select Case [value]

Case "1"
  area_type = "aapent omraade"
Case "11"
  area_type = "bebygd"
Case "12"
  area_type = "samferdsel"
Case "21"
  area_type = "fulldyrka jord"
Case "22"
  area_type = "overflatedyrka jord"
Case "23"
  area_type = "innmarksbeite"
Case "30"
  area_type = "skog"
Case "50"
  area_type = "aapen fastmark"
Case "60"
  area_type = "myr"
Case "70"
 area_type = "isbre"
Case "80"
 area_type = "vann"
Case "99"
 area_type = "ruralt (ikke kartlagt)"
End Select


texts:
When I run ./configure --with-grass=/usr/local/grass-6.5.svn/ I get the following errors:

..........
checking for gdal-config... /usr/bin/gdal-config
using /usr/lib/gdalplugins as GDAL shared library autoload directory
checking for G_asprintf in -lgrass_gis... no
configure: error: --with-grass=/usr/local/grass-6.5.svn/ requested, but libraries not found!  Perhaps you need to set LD_LIBRARY_PATH to include /usr/local/grass-6.5.svn//lib?

Also if I use --with-grass=/usr/local/grass-6.5.svn, I get the following error:

checking for gdal-config... /usr/bin/gdal-config
using /usr/lib/gdalplugins as GDAL shared library autoload directory
checking for G_asprintf in -lgrass_gis... no
configure: error: --with-grass=/usr/local/grass-6.5.svn requested, but libraries not found!  Perhaps you need to set LD_LIBRARY_PATH to include /usr/local/grass-6.5.svn/lib?

I have created a grass.conf as instructed elsewhere inside /etc/ld.so.conf.d/ and have added following lines in it:

/usr/local/lib
/usr/local/grass-6.5.svn/lib

I ran sudo ldconfig after that too.

I came across a similar issue posted in the GRASS-Dev mailing list. As an answer it was instructed to: You could just generate the symlinks yourself, e.g.:

    cd /usr/local/grass-6.5.svn/lib
    for file in libgrass_*.so ; do ln -s $file ${file%.6.5.svn.so}.so ; done html.

I ran,

cd /usr/local/grass-6.5.svn/lib
ln -s $file ${file%.6.5.svn.so}.so

It made no difference to the errors. Did I make any mistake? What do I need to do?


texts: Ogr Python API: set null values in feature

texts:
Im trying to convert OGR geodata from one driver to another, like the command line tool ogr2ogr do, but using the python API.
My problem is, that the NULL values in the input data are converted as "None" string in the output data.

The Ogr SetField method:

output_feature.SetField('linetype',None)

does not set the attribute 'linetype' =  None (or null) in the output feature, as expected, but "None" as string.
I tested this with PostgreSQL and GeoJson driver, but probably affected all drivers.

Complete code sample from here:

outDriver = ogr.GetDriverByName("GeoJson")
outDataSource = outDriver.CreateDataSource('test.geojson')
outLayer = outDataSource.CreateLayer("test_layer")

# Add an name field
idField = ogr.FieldDefn("name", ogr.OFTString)
outLayer.CreateField(idField)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(1198054.34, 648493.09)

# Create the feature and set values
featureDefn = outLayer.GetLayerDefn()

#feature 1
feature_1 = ogr.Feature(featureDefn)
feature_1.SetGeometry(point)
feature_1.SetField("name", 'name_1')
outLayer.CreateFeature(feature_1)

#feature 2, name attribute should be NULL but is saved as "None" (string)
feature_2 = ogr.Feature(featureDefn)
feature_2.SetGeometry(point)
feature_2.SetField("name", None)
outLayer.CreateFeature(feature_2)

# Close DataSource
outDataSource.Destroy()


texts: QGIS-bin has stopped working; QGIS crashes at startup

texts:
How can I recover a functional version of QGIS, and hopefully also recover the project I'm working on? I'm using Windows 7 64 bit. Here is the screenshot.

I have read the threads to other similar questions and haven't resolved the problem.

A new QGIS installation, removing 2.18 La Palmas and replacing with 2.18.11, also did not work
Removing C:\Program Files\QGIS 2.18\apps\qgis\processing folder did not work

This crash happened as I was trying to add svg markers to my composer - as a final touch on my completed project! There are just no north arrows (but I did find them in QGIS here: C:\Program Files\QGIS 2.18\apps\qgis\svg\symbol when I tried to add some svg markers I had downloaded. They just don't show up under the svg available markers for arrows in composer? Anyway, something I did in that effort crashed the whole program and all future installations.


texts: Connecting to Oracle User Schema Geodatabase using ArcObjects?

texts:
I've set up an Oracle user Schema geodatabase and loaded some data into it. I now want to connect to it programmatically using ArcObjects, the problem is I don't see how. Where can I set the newly created schema name in my connection property set, or how do I change the schema value after I have connected? In ArcCatalog it is done in the Geodatabase connection properties window :

The code below works fine connecting to the default SDE schema. What do I have to change to connect to a user geodatabase? I have tried examining IDatabaseConnectionInfo2 and IWorkspaceProperties but I don't see where I can set or change the schema that is connected to.

 public static IWorkspace ConnectToTransactionalVersion(String dbclient, String dbConnProp, String user, String password, String database, String version, String authentication)
    {
        IPropertySet propertySet = new PropertySet();
        propertySet.SetProperty("DBCLIENT", dbclient);
        propertySet.SetProperty("DB_CONNECTION_PROPERTIES", dbConnProp);
        propertySet.SetProperty("DATABASE", database);
        propertySet.SetProperty("USER", user);
        propertySet.SetProperty("PASSWORD", password);
        propertySet.SetProperty("VERSION", version);
        propertySet.SetProperty("AUTHENTICATION_MODE", authentication);
        Type factoryType = Type.GetTypeFromProgID(
            "esriDataSourcesGDB.SdeWorkspaceFactory");
        IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance
            (factoryType);
        return workspaceFactory.Open(propertySet, 0);
    }

I tried adding "SCHEMA" as a property in propertySet and it gets ignored, I also tried setting the version property to the new user schema default version but it errors with "version not found"


texts: Using QGIS vs. SVG to export with area names?

texts:
Is there the possibility of exporting a SVG files (from shapefile) which retains the area names?

Basically any closed path should have the original name attributed to shapefiles. But not only. I would also have individual geoid points appointed in the same way. Even so unique, "point_AREANAME", while the area should be called simply "AREANAME".

Is this feasible?


texts: 'ObjectId field is missing' error when using Esri Leaflet updateFeature

texts:
I am trying to update a feature on the map upon click of a button. I am using updateFeature function of esri leaflet feature layer.

Here's the code that I have written:

layer_view_18.updateFeature({
    type: 'Feature',
    id: target.feature.properties.ID,
    geometry: geometry,
    properties: markerProperties
}, (err, response) => {
    console.log('response from updateFeature -- ', err);
    console.log('response from updateFeature -- ', response);
    // layer_view_18.redraw(target.feature.properties.ID)
});

The console show error saying 'ObjectId field is missing' and the response is undefined.

The map service layer has OBJECTID field defined in it along with other fields.

Has anyone faced this issue before?

I am following the editing example from esri-leaflet documentation.

Is there something that needs to be changed in the code?


texts: Getting intersecting polygons of points using OGR and Python

texts:
I have multiple files that essentially have anywhere from 2000-5000 points each.
I have a polygon file that has about 10000 polygon shapes.
For each file containing the points, I want to get all the polygons that the file intersected.
Checking each polygon for each point and using Contains() takes too long. 2000 x 10000 iterations for each file essentially.
My current approach is to create a polygon around the points then get the intersections of that one bounding polygon to the 10000. This is fast but leads to me picking up intersecting polygons that don't actually contain the points but intersect the bounding polygon I had made around the points. So I have to go back then and do another check for intersections that don't actually contain points.
#go through files containing points and create a bounding polygon around each
for file in point_files:
        pts_ds = ogr.Open(file)
        pts_lyr = pts_ds.GetLayer()
        geom_collection = ogr.Geometry(ogr.wkbGeometryCollection)
        for feat in pts_lyr:
            geom_collection.AddGeometry(feature.GetGeometryRef())
        bounding_polygon = geom_collection.ConvexHull()
        #then save bounding polygon to file

#open up each of the bounding polygons I had made and intersect it against the comparison polygons
comparison_polygons_ds = ogr.Open(comparison_polygons.shp)
comparison_polygons_lyr = comparison_polygons.GetLayer()
matches = {}
for file in bounding_polygon_files:
    b_polygon_ds = ogr.Open(file)
    b_polygon_lyr = b_polygon_ds.GetLayer()
    for b_polygon_feat in b_polygon_lyr:
        b_polygon_geom = b_polygon_feat.GetGeometryRef()
        for comparison_polygons_feat in comparison_polygons_lyr:
            comparison_polygon_geom = comparison_polygons_feat.GetGeometryRef()
            if b_polygon_geom.Intersects(comparison_polygon_geom):
                matches.setdefault(comparison_polygon_geom.GetField(“ID”), []).append(file)

The end result I want is essentially: file1 intersects x,y,z polygons. file2 intersects a,b,c polygons.
How do I achieve that?
I am unfamiliar with GIS and using Python with the OGR package.


texts:
I have a WMTS service which I want to use for my work.

The problem is that I can't use it with Autocad Map because it consumes only WMS, or ArcGIS 10, because ArcGIS 10.1 and above uses WMTS.

Is there any possibilities to convert it on the fly in a WMS?


texts: How to show line symbol (arrow) independently of line length in ArcGIS?

texts:
I am manually creating a polyline layer in ESRI geodatabase. Lines are of 2 vertices (start-end), variable length but all very short. I am using directional symbol with arrow at the end (or middle). And I have a small annoying problem: I cannot see the shortest lines. The arrow just disappears when its size is less than the line at that scale. What I get is this:

And I have troubles to see short lines like the upper two especially when I add more layers to the map, because I am creating the lines based on other geodatabase layers.

Increasing the width is not good because it shows a perpendicular expansion of a real line. I just want ArcGIS to show me the arrow symbol always, no matter the visible length of the line. I was searching in Symbol properties as it was advised also here
Showing directional line symbol in ArcGIS for Desktop?
but could not find the solution. When I google, I find only opposite situations (people want to show different symbols depending on the length). Converting line ends to points to show symbols is not possible because I am editing the layer and adding new lines. So I really need to see them well, and their direction. What I want is this:

I am using ArcGIS 10.2 for Desktop. Maybe this is a fault of ArcGIS? Because QGIS is showing arrows as I need:


texts:
We cannot seem to find any way to export our projects as PNG/JPEG through the print layout/composer. We draw normal line and point SHP files over the Bing/Google Earth etc layers, but the exported image is simply blank with only our line and point layers showing.

We can sometimes get exports succesfully when the scale is over 3000, but due to the detailed information we need to depict, we need to get down to a scale of 1000.

Does anyone have any solutions or alternative options we can try?


texts: LIBKML using Python very slow - using OGR to transform Shapefile to KMZ/KML

texts:
I am trying to convert a shapefile to a KML (or KMZ) using Python, and OGR and the LIBKML driver. The process appears to hang, with the kml/z created, but sitting at 1kb for 5+ minutes.
The input shapefile has about 20 columns that are a mix of Double, and Text, and is about 8mb in size.
Any bright ideas?
The following code fails/hangs:
in_ds = ogr.GetDriverByName('ESRI Shapefile').Open(os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.shp')
out_ds  = ogr.GetDriverByName('LIBKML').CopyDataSource(in_ds, os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.kmz')

The following code works, and completes in about 1.5 seconds.
in_ds = ogr.GetDriverByName('ESRI Shapefile').Open(os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.shp')
out_ds  = ogr.GetDriverByName('KML').CopyDataSource(in_ds, os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.kml')

Worked Example:
from osgeo import gdal
from osgeo.gdalconst import *
from osgeo import osr
from numpy import *
from osgeo import ogr

out_dir = 'C:\Temp\'
output_dir = 'C:\Temp\'

print('Writing Shapefile')

in_ds = ogr.GetDriverByName('ESRI Shapefile').Open(os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.shp')
out_ds  = ogr.GetDriverByName('LIBKML').CopyDataSource(in_ds, os.sep.join([output_dir, os.path.basename(out_dir)]) + 'impact.kmz')


texts: Converting KML to layer in ArcGIS 10 - nothing happens

texts: I have converted many KML files to layer in ArcGIS many times before without any problem.  I just did a conversion (Arc Toolbox>Conversion tools> KML to layer) and the file is showing up in the table of contents but there is no layer. Any ideas or suggestions??


texts: Modify attributes of WFS layer in OpenLayers 3

texts:
I need to modofy attributes of features appeared on a map using OL3.
The features live in a postigs db and they are serverd through Geoserver.
What I want is to modify some of the non spatial attributes of each feature.

What I did in order to accomplish this is to create a menu in which the user selects the values he wants to modify for a specific feature and then save the changes with a button. On button click I would send an ajax request (jquery/php) in which I would update the table values based on the fid of the table.

If I understand correctly using WFS-T is a more appropriate way. Is that right? Can I change non spatial attribute values with WFS-T? If yes, is there an example using OL3?


texts: QGIS, MMQGIS, join by attribute function, cannot read the text(tab delimited file)

texts:
I am trying to join none spatial data with geographic data using the MMQGIS tool, on QGIS 3.14. However, I am running into an issue, when I try to carry out the function the MMQGIS tab will not let me select any fields for the Text(tab delimited) file. As such, I cannot carry out the function, due to the fact that I cannot select a field to match. It is worth noting when I import the text(tab delimited) file as a layer it works fine, I can see all the fields in the attributes table and all the data is there. Also note this is part of a course the the data set has been used by 100s of students, so it is unlikely that the data set is the issue.
I have tried using a CSV file in place of a text (tab delimited) file but to no success. It does let me select fields but only fields 1,2,3,4 and so on not the actual names of the fields, and the operation still fails. Please refer to the attached photo for further details, this should be easy as I am following a step by step tutorial and it still does not want to play ball.
I have tried it on to different computers one running QGIS v 3.10 and I have tried it on Qgis v2.8.
What should I do from here?


texts: How to use multiprocessing or pool to boost up the performance of ArcPy script

texts:
I created some arcpy script which actually runs slow on huge dataset. such as 100000 feature in it.
I read from below path that multiprocessing can boost up arcpy script script performance.
 http://blogs.esri.com/esri/arcgis/2011/08/29/multiprocessing/
I tried it myself and applied it in my code.

It is giving an error "Could not find file: from multiprocessing.forking import main;main().mxd"

I could not understand what does this error means,
and was unable to get any reference or help from any Internet source about the error.

Below is the code which actually find duplicates in a neaby buffer area(defined by user input)

import arcpy
import os
import multiprocessing
# defining input
inputFeatureClass = arcpy.GetParameterAsText(0)
InputField = arcpy.GetParameterAsText(1)
Distance = arcpy.GetParameterAsText(2)
#checking whether UNIID  and dup field exist or not, if not than adding it
to_add = ["UNIID","Dup"]
#Create a list of existing field names
fieldList = arcpy.ListFields(inputFeatureClass)
fieldName = [f.name for f in fieldList]
for field in to_add:
  if not field in fieldName:
    arcpy.AddField_management(inputFeatureClass,"UNIID","TEXT")
    arcpy.AddField_management(inputFeatureClass,"Dup","SHORT")
#updating squence number into UNIID Field
rows = arcpy.UpdateCursor(inputFeatureClass)
counter = 1
for row in rows:
    row.UNIID = counter
    rows.updateRow(row)
    counter += 1

#setting up progressor
RowCount  = arcpy.GetCount_management(inputFeatureClass)
PACnt = int(RowCount.getOutput(0))
arcpy.SetProgressor("step", "Finding Duplicate...",0, PACnt, 1)

def Find_dup(oid_range):
  i = oid_range[0]
  j = oid_range[1]
  lyr = arcpy.management.MakeFeatureLayer(inputFeatureClass, 'layer{0}'.format(i),"OID >= {0} AND OID <= {1}".format(i, j))
#selecting neaby Point for every point inside searchcursor
  with arcpy.da.SearchCursor(lyr, ["UNIID", InputField]) as cursor:
    for row in cursor:
        value=row[0]

        field="UNIID"
        exp=field + "='" +value+"'"
        point_lyr=arcpy.SelectLayerByAttribute_management(lyr,"NEW_SELECTION",exp)
        arcpy.SelectLayerByLocation_management(point_lyr, overlap_type="WITHIN_A_DISTANCE", select_features=point_lyr, search_distance=Distance + " Meters", selection_type="NEW_SELECTION")
#pushing selected points into array
        occursTwice = []
        cursor = arcpy.da.SearchCursor (lyr, [InputField])
        for row in cursor:
            occursTwice.append (row[0])
        del cursor
#counting occurance of every Field Value inside Array and updating into dup field
        cursor = arcpy.da.UpdateCursor (lyr, [InputField, 'Dup'])
        for row in cursor:
# check whether dup contains already a duplicate value, if yes than do not update lower value
            if occursTwice.count(row[0]) > row[1]:
              row[1] =  occursTwice.count(row[0])
            cursor.updateRow(row)
        del cursor
        del point_lyr
        arcpy.SetProgressorPosition()

def main():
  oid_ranges = [[0, 1000],[1001, 2000], [2001, 3000], [3001, 4000], [4001, 5000],[5001, 6000], [6001, 7000], [7001, 8000], [8001, 9000], [9001, 10001]]
  arcpy.env.overwriteOutput = True
  pool = multiprocessing.Pool(processes=4)
  pool.map(Find_dup, oid_ranges)
  pool.close()
  pool.join()
if __name__ == '__main__':
  main()

#Deleting UNIID Field and Unselect everything
arcpy.DeleteField_management(inputFeatureClass,"UNIID")
arcpy.SelectLayerByAttribute_management(inputFeatureClass,"CLEAR_SELECTION")
arcpy.ResetProgressor()


texts: Writing a raster file from multi-dimension array using Python

texts:
I've been attempting to create a raster file that has multiple bands (10 or more) from a numpy array with multiple dimension(i.e. 10 or more). The script I'm currently using exports the TIFF file I want to produce but the file turns out to be black (blank) image.

The shape of the data has a dimension of (1120,872,12) (rows, cols, bands) and I wanted to create a raster file that has 12 bands.

input_array.shape
 (1120, 872, 12)

def CreateTiff(output, array, driver, noData, GeoT, Proj, DataType):

    rows = array.shape[0]
    cols = array.shape[1]
    band = array.shape[2]

    array[np.isnan(array)] = noData
    driver = gdal.GetDriverByName('GTiff')

    DataSet = driver.Create(output, rows, cols, band, gdal.GDT_Float32)
    DataSet.SetGeoTransform(GeoT)
    DataSet.SetProjection(Proj)

    for i, image in range(array.shape[2], 1):
        DataSet.GetRasterBand(i).WriteArray(image)
        DataSet.GetRasterBand(i).SetNoDataValue(NoData)
    DataSet.FlushCache()
    return Name

with this function,

# the raster output layer

output_file = r'rasterfile.tiff'

# Writes raster

CreateTiff(output_file, input_array, driver, noData, GeoTransform, Projection, DataType)

As I mentioned above, this function creates the raster file I was hoping to get but when I inspect it in QGIS, the file is a black(blank) image and the values are zeros.

Any thoughts?


texts: Get value from applied definition query and update field value

texts:
I have a layer "ABC", in which I have a field [LABEL]. In the field [LABEL]  some feature have been filled with the value "TSL/VNR/DT8A"
    definition Query is already set to [LABEL] ='TSL/VNR/DT8A'

I now need a script that could read the content of the definition query, pull out the single quoted text and fill this value to other feature.

Below is the code I am preparing.

mxd = arcpy.mapping.MapDocument("CURRENT")
layers = arcpy.mapping.ListLayers(mxd,"ABC")
# Getting Current Definition Query
DEFQUERY =  layers[0].definitionQuery

# Print DEQUERY
# [LABEL] ='TSL/VNR/DT8A'

***# challange is to extract single quoted text from entire definition query.***
SingleQuotedText = ***some expression to pull out single quoted text from DEFQUERY***

#releasing definition query.
lyrnm = arcpy.mapping.Layer("ABC")
lyrnm.definitionQuery = None

#cursor to update unique value in field LABEL
cursor = arcpy.UpdateCursor("ABC")
for row in cursor:
        row.LABEL = SingleQuotedText
        cursor.updateRow(row)

My challenge is to pull out the single quoted text and put it into the cursor.
I am not having good knowledge of regular expression. It may help me, but don't know how to apply it.


texts: How do I get the value of the pixel under the cursor?

texts: Is there a tool or something that will give me the value of the pixel under the cursor? What I've got in mind in something like the eyedropper tool in graphics programs that allows you to interactively sample points in a raster image.


texts:
I am trying to delineate the riparian corridors using GRASS in R. But I am having some issues.

So is there a plugin in QGIS using which I can extract riparian corridors using the already available upstream corridors?


texts: How exactly does QGIS compute the 3D Tiles for rendering a 3D Map View

texts:
In QGIS, we can click "View > New 3D Map View". This allows us to see our map in 3D, as long as we have a DEM layer available.
But QGIS takes some time to actually compute the 3D tiles here:     ↓

Before they render:

Here, I went "under" the ground with the camera, it's a trick to show a little of the tiled structure of the terrain. Colors represent the slope, the steeper, the more red it is.
I would like to figure out the underlying technical details on how QGIS does compute these 3D tiles, what is the algorithm involved, are these tiles computed in, e.g. the Cesium 3D Tiles format, or is this an internal QGIS format? There is not much information during this computation, that's why I'm asking. And I'm not able to get any result on the QGIS documentation itself at: https://docs.qgis.org/3.16/en/search.html?q=3d


texts:
I'm trying to setup Geonode and Geoserver advance security. I'm currently doing setup for the Geoserver.

I'm stuck at this step Setup of the GeoNode REST Role Service.
I cannot add new role service as i did not have AuthKEY REST - Role service from REST endpoint in the role service option.

AuthKEY REST - Role service from REST endpoint.

currently using

Django (1.6.11)
Geoserver 2.4.7
GeoNode (2.4.dev20170302183855, /vagrant/geonode)
django-oauth-toolkit (0.7.2)


texts: Getting distance from point to nearest line string

texts:
I am trying to use st_distance to find the nearest line string to a point. I have two tables in cartodb one is a table of points and the other table is poly lines. To see if I had the logic right I did a st_distance query on the points table from a lat,long point. Everything worked as planned. Then I tried creating a query that would find the closest polygon from a point in my point table.

The error I am getting is this "cross-database references are not implemented: cartodbpoint.the_geom.stdistance" and the query I am trying out is:

    SELECT cartodbpoint.cartodb_id,
    MIN(cartodbpoint.the_geom.STDistance(cartodbline.the_geom)) as MinDistance
    FROM cartodbpoint, cartodbline
    GROUP BY cartodbpoint.cartodb_id

My goal is to get the distance from that point to the nearest line string.


texts: I was wondering if there is a way to convert heights that are logged from a normal bluetooth GPS to AHD heights.


texts:
I am trying to create a geoTIFF from an .hdf MOD35_L2, I am trying either a Matlab code, either I heard something about LDOPE-1.7 + MRTSwath but I am not really succeeding.

Have anyone ever tried to do so?


texts: Regrouping points together based on a sum

texts:
I have a vector layer containing points. Each point has an attribute called 'number of houses'. I want to dissolve this vector layer in many areas, each area must contain between 643 and 653 house and points must be near to each other.

The result must be something like this:

I work with QGIS and Python.

I am new in QGIS and Python and I don't know how to enter my shapefile or the coordinates as an input in those algorithms. And I don't want to divide my set into equal clusters, but into clusters based on the numbers of houses which are written near the point.


texts:
Lacking much programming experience, I've spent a lot of time familiarizing myself with HTML, CSS, and SQL over the past couple of days so that I could make a layer selector for a CartoDB visualization.

I eventually made a layer selector, but my SQL queries weren't working for some reason. Apparently I had to make the data table non-private. I was informed that I could have looked at the console and it would have told me the error.

Can someone tell how I can access this console?

I have no idea what that is. I was thinking it could be accessed via Linux terminal, but my platform right now is Windows.


texts:
ArcGIS supports building polygons within the flex component that contain self-intersecting rings, as well as rings that touch at multiple points (breaking the "arbitrary" rules here: http://postgis.net/workshops/postgis-intro/validity.html ).

This causes incompatibility when passing these shapes to PostGIS.  The article above mentions buffering to fix some of these issues, but this seems to only fix some cases, and often ends up with a different shape polygon than was intended (eg, only half of the poly).

Is there a way to convert these to monotonic polygons within flex so that this problem is eliminated?  Or a Java API to do the same?

Or is there some other approach entirely to solving this problem?


texts: Finding distance to nearest point from feature using Google Earth Engine?

texts:
In Google Earth Engine I have a FeatureCollection of 14,002 points. For each of the points I want to find the next nearest point's distance away. I don't need to record the nearest point ID, only the distance to it and save it in the feature. I haven't been able to find anything in the Docs that would do this.

Essentially I am trying to implement ArcMap's NEAR tool (found here: http://desktop.arcgis.com/en/arcmap/10.3/tools/analysis-toolbox/near.htm
into something similar in GEE but haven't found a way yet.


texts: Concatenating list of strings using ArcPy?


texts: Where are the temporary output layers from QGIS processing algorithms stored?

texts:
Context

I have always worked with memory outputs from several QGIS geoalgorithms within the Processing Toolbox. If it wasn't specified a location when saving the result, the output was stored as a temporary layer in a specified path, which was in general something like:

  C:\Users\xxxxx\AppData\Local\Temp\processinge252173536ab4ef987fcb13829939857\c35bd5c897f24c5eb13b9aad985b9652\OUTPUT.shp

If the temporary layer wasn't saved before quitting QGIS, it was automatically deleted from the disk.

Issue

Using QGIS 2.18, I casually noted that I'm no longer able to find the temporary output layers in the Temp folder: this behavior appears to only belong to the QGIS geoalgorithms (instead, the temporary outputs from SAGA or GRASS algorithms are still saved in the Temp folder, as usual).

If I select one output from the QGIS geoalgorithms and go to Properties > General, I see something like this in the Layer source field:

  Polygon?uuid=29cdbb21-d299-4b01-b735-0004971e1203&crs=EPSG:32632&field=id:integer&field=value:string&field=area:double&field=perim:double

As you can see, there isn't any reference to a path in my disk, so I can't find the physical position of the output.

In Settings > System, I have this path for the temporary files:

C:\Users\xxxxx\AppData\Local\Temp

In the Processing options, instead, I have:

C:\Users\xxxxx\.qgis2\processing\outputs

There is something that I don't know on how the temporary output layers are managed since the current version of QGIS? I encountered this issue using QGIS 2.18, but I don't exclude that it could be started from previous versions of the software.


texts: JTS: generate a polygon geometry from line geomtry

texts:
I have a line and I wish to convert it to a polygon.
The point is this line is NOT a ring
it's an straight Line. And I want to use an specific width to extend line horizontally and create a polygon.

Using JTS.
is it possible? How can I make it happen?


texts: Adding a Tool to a DockableWindow add-in

texts:
I have inherited a C# ESRI (ArcMap) add-in that needs a UI tweak (the rest of the changes are completed).  And I had no previous experience with such add-ins before starting on this.

The add-in in question has a DockableWindow, but I need to be able to click on the map and get the coordinates of the click (the inherited version required typing in the coordinates manually).

Research indicated that the only way to do this was to create a new class inheriting from a Tool, and use OnMouseDown()/OnMouseUp() to do the work.  (I did try overriding OnMouseClick() in the DockableWindow, but that seems to only capture clicks within the DockableWindow.  Not a huge surprise there.)

In any event, I created the Tool, added a line for it in the .esriaddinx file, and saw that it was added to ThisAddIn.IDs.

Everything I've read seems to indicate that the next steps are to get the UID for it, do a Find on the document's command bars using that UID, then set the app's currentTool to my tool.

The only problem is that, when I try this, the Find always returns null.

I've even tried creating my own command bar, adding the Tool to that command bar, displaying the command bar, and doing the Find on that command bar.  Still, I always get null.

Is there another step that I'm missing?  (And sorry about the lack of actual code; the development machine is not near a computer connected to the internet.)


texts: Read/Write generic tiff tags using GDAL API

texts:
Is it possible to read/write generic tiff tags with GDAL?  I know some tags are supported (http://www.gdal.org/frmt_gtiff.html), but is there a way to access other tags?

For example, there are geo tags I would like to access directly.  If this is not currently supported are there future plans?


texts: How to add a new data frame to an MXD using C#

texts:
Does anyone have code in C# to add a new data frame (map) to an MXD?  I have an MXD created in ArcMap 10.1 that has an existing data frame.  I need to programmatically add a new data frame and populate it with several layers from a personal geodatabase.  I'm using VS2010 with .Net 3.5.  I have created a MultiItem that reads the names of folders that contain the PGDB files.  After selecting a name from the list (which is populating correctly) it triggers OnSelChange.  In that routine I need to create a data frame named based on the folder name selected from the MultiItem list.  The name comes in just fine, but I'm stuck creating the new dataframe.  I tried using IMaps.Create() but it doesn't add anything to the count of maps.  I have tried then using IMaps.Add() but VS reports that Add is not defined for IMaps contrary to being described in the ArcObjects reference.

The only procedure I have found is the one written in VBA using ArcObjects 10.  It seems rather convoluted given there are supposed to be functional methods for IMaps to do exactly what I need.  Any help would be greatly appreciated.


texts:
I am starting to research on the topic related to use of GIS in public health. I have started with the Wikipedia article.

I would like to apart from epidemic spread of disease, what are the other application areas? Currently active research groups/universities working on these topics? World Health Organization projects? Resources that would help students researching in this area?


texts:
I want to do the following:

from a point layer's table,

1) grab all GEBID (an attribute) and create a list(1) with unique GEBIDs

2) iterate through the list(1) in a loop:

  2a) stuff all features with matching GEBIDs into a temporary list(2)

  2b) pick the entry with the largest Area (another attribute)

  2c) write the entry into a third list(3)

3) use the third list(3) to select the features

The resulting selection should be like the red text:

Here's my code so far, but I'm stumped at 2b)..

import arcpy

IDList = []
rows = arcpy.SearchCursor("A_Gebauede_Dachform_F2P")

for row in rows:
    IDList.append(row.getValue("GEBID"))

IDSet = set(IDList)
IDSetList = list(IDSet)

for i in range(len(IDSetList)):
    rows = arcpy.SearchCursor("A_Gebauede_Dachform_F2P","GEBID = '"+IDSetList[i]+"'")
    ???


texts: 1:250 streams data - USA wide - where can I find it?

texts:
I have 1:250 streams data for a few western states. I find this scale to be excellent in many settings as base hydro data.

I have never been able to find a USA or better yet, NA wide dataset of this scale. Does anything like this exist for the USA or NA?

Look at the attached Idaho jpg for an example of this data


texts: Can I host multiple map services from a single file geodatabase?

texts:
Can I host and consume multiple map services from a single file geodatabase?

The data will be view only. Feature service capability will not be turned on, so no editing will be performed.


texts: How to call function inside another function in qgis python plugin?

texts:
I am developing qgis plugin.

class HyperspectralClassifier:
 """QGIS Plugin Implementation."""

  def __init__(self, iface):
  .
  .

    self.dlg.pb8.clicked.connect(self.plotspect)
  def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError:
        raise ValueError("window_size and order have to be of type int")
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    # precompute coefficients
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    # pad the signal at the extremes with
    # values taken from the signal itself
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
    return np.convolve( m[::-1], y, mode='valid')

  def plotspect(self):
   #QMessageBox.information(None, "Clicked coords", " Row: " + str(row) + " Column: " + str(column) )
   iinf=self.dlg.iin.text()
   img=gdal.Open(iinf)
   inputArray=img.ReadAsArray()
   ia1=inputArray[4:195,141,194]
   ia1=savitzky_golay(ia1, 51, 3)
   plt.plot(ia1,label="Original Spectrum")
   plt.show()

but it  gives error like this: global name 'savitzky_golay' is not defined
what is the problem?


texts: Making ring buffers (donuts) in Google Earth Engine

texts:
I am trying to generate multiple buffers around a layer that represents cities (FeatureCollection). Use the following code to create the buffers:
// Function to create buffers of different sizes in Google Earth Engine .
The problem is that this code generates polygons that include the inside area. I need to generate buffers on the outside, like doughnuts or perimeter rings.
I tried to use .difference to remove the original polygons from the buffers but I get the error "....difference is not a function"
Can anyone help me?
var geometry =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      },
      {
        "type": "rectangle"
      },
      {
        "type": "rectangle"
      }
    ] */
    ee.FeatureCollection(
        [ee.Feature(
            ee.Geometry.Polygon(
                [[[-92.93044076036253, 43.38666530027385],
                  [-92.93044076036253, 43.371630789523294],
                  [-92.90812478135862, 43.371630789523294],
                  [-92.90812478135862, 43.38666530027385]]], null, false),
            {
              "system:index": "0"
            }),
        ee.Feature(
            ee.Geometry.Polygon(
                [[[-93.016614771593, 43.389328699363205],
                  [-93.016614771593, 43.36986507859759],
                  [-92.99155221055784, 43.36986507859759],
                  [-92.99155221055784, 43.389328699363205]]], null, false),
            {
              "system:index": "1"
            }),
        ee.Feature(
            ee.Geometry.Polygon(
                [[[-93.09137330125608, 43.39935439198635],
                  [-93.09137330125608, 43.38213974024519],
                  [-93.06734070848265, 43.38213974024519],
                  [-93.06734070848265, 43.39935439198635]]], null, false),
            {
              "system:index": "2"
            })]);

var RingBuffer = function (Poligono, DistBuff) {

var BufferFeature = function(f) {
  f = ee.Feature(f);
  var buffer_size = f.get('buffer_size');
  return f.buffer(buffer_size);
};

var BufferFeaturesByDistance = function (fc, buffer_size) {

var SetBufferSize = function(f) {
    return f.set({'buffer_size': buffer_size});
  };

return (Poligono).map(SetBufferSize).map(BufferFeature);
};

var Buffer = BufferFeaturesByDistance(Poligono, DistBuff)

return (Buffer).difference(Poligono)
};

var Buff1 = RingBuffer (geometry, 2000);

Map.addLayer(Buff1);


texts: Zooming & panning to selected objects using ArcGIS Engine?

texts:
I'm writing a program using ArcGIS Engine, and need to zoom & pan into the map to show the selected content.

Loading and displaying the map does work using something like this:

AxMapControl _mapControl;

// in constructor:
_mapControl = new AxMapControl();

// in loading
_mapControl.LoadMxFile(@"C:\Users\me\Documents\TestProject.mxd");

Afterwards I select some objects and want to zoom into them, but I can't find a way to get the selected extent - any pointers to achieve this?


texts: Snapping of ArcGIS Online features fails in ArcGIS pro

texts: I'd like to edit a feature-layer that is published in ArcGIS Online with ArcGIS pro. I urgently require the snapping function. However, after enabling snapping I can only snap to the feature that is edited but not to other features of the layer. Does anyone has an idea how to solve this problem?


texts: I had a point file of an area around 70 km^2 collected through field survey using Auto Level. The survey was carried out in the lines 25 m apart. I converted the readings from auto level into point shapefile and then converted the point shapefile into DTM using Grid- IDW (nearest neighbour selection) tool with search radius of 25 m in QGIS. However, in the DTM prepared through GIS has lot of discontinuity in form of strips at the distance of 25 m. I wish to smoothen these strips. Can you suggest the reasons?


texts: Making ArcGIS Attachments compulsory (required) for each feature added?

texts: Is there a way to make attachments compulsory. That is, each time a user adds a new point or record, they have to include an attachment. I am exploring the Relationship class (Changing rules within it) but it looks like it can only be applied to subtype. There is no subtype in my feature class. I want to make sure each time a record comes in, there is a picture tag to it.


texts: How to project WRF data (lat/lon) onto ArcGIS maps

texts:
I have climate data generated using the WRF model, in lat/lon format (lambert conformal projection), in an Excel file.

I'm looking to project / overlay this data onto a map of the UK using ArcGIS, however the projection is only showing as a single point, and not over the UK!

Has anyone got any experience of overlaying WRF output (in lat/lon format) onto a GIS map?


texts: QGIS : cannot change the z-coordinate of a point shapefile using the wkb geometry format

texts:
I am using QGIS Desktop 2.5.0-Master with GDAL/OGR 1.11.0 on a windows 7.

I am working on a little pyQGis script that should modify the z-coordinate of points in a vector shapefile based on a value read in a raster (DEM). I really want to edit the geometry and not to create/update a feature attribute that hold the z-coordinate.

When executing the code below, the execution of the script goes through without throwing any error, but it seems that the changes to the wkb geometry are not committed...

I went then into the QGIS python console to directly check the outputs of my for loop on the selected features in the vector layer. What I saw is that the call to "geom_ogr.SetPoint(0,newX,newY,newZ)" was really resulting in a change of the geom_ogr.GetX(), geom_ogr.GetY() and geom_ogr.GetZ() values. However, once out of the for loop, I can't commit those changes in the shapefile...

Is there a reason for not being able to write permanently new geometry information in a point vector layer ??

--------PLUGIN CODE---------

vector_layer=vector ->definition of a parameter used in the GUI proposed by QGIS for tools
raster_layer=raster ->definition of a parameter used in the GUI proposed by QGIS for tools

from qgis.core import *
from PyQt4.QtCore import *
from osgeo import ogr

vector_path = vector_layer
raster_path = raster_layer

vector = QgsVectorLayer(vector_path,"vector_lyr","ogr")
if vector.isValid()!=True:
    print("impossible de charger la couche vecteur")

raster = QgsRasterLayer(raster_path,"raster_lyr")
if raster.isValid()!=True:
    print("impossible de charger la couche raster")

selection = QgsFeatureRequest()
selection.setFilterRect(raster.extent())

vector.startEditing()

for f in vector.getFeatures(selection):
    wkb = f.geometry().asWkb()
    geom_ogr = ogr.CreateGeometryFromWkb(wkb)
    altitudeDic raster.dataProvider().identify(QgsPoint(f.geometry().asPoint().x(),f.geometry().asPoint().y()),QgsRaster.IdentifyFormatValue)
    z = altitudeDic.results()[1]
    geom_ogr.SetPoint(0,geom_ogr.GetX(),geom_ogr.GetY(),z)

 vector.commitChanges()


texts: i.atcorr [GRASS] - has anyone used it and actually had success?

texts:
I am trying to apply i.atcorr to a landsat mosaic with GRASS. These are my 6S parameters:

18
01 31 22.10 178.940 -17.263
1
1
0
0.115
-0.270
-1000
 116

Everytime I run it, it returns me a raster with only one value. I've changed the values, I've tried with dummy data, I've read the documentation several times. I have read literally all the questions here in stack gis about the use of i.atcorr with GRASS, and in their mailing list as well. There seems to be a lot of people that has the same output, this sort of blank raster, but there is not a clear solution to this.

Has anybody applied this algorithm successfully to Landsat 8 and, if so, any idea on I am possibility doing wrong?

**Edit:
GRASS GIS 7.4.4. on Windows 10.


texts: Accessing fields of created features in your PyQGIS

texts:
I'm new using PyQGIS 3 and I'm wondering how to access the field information from a layer generated by script.
Inside my processAlgorithm function I have the following:
alg_params = {
    'EXPRESSION': parameters[self.IN_EXPRESSION],
    'INPUT': parameters[self.IN_SHP],
    'OUTPUT': parameters[self.OUT_SHP]
}

final_output = processing.run('native:extractbyexpression', alg_params, context=context, feedback=feedback, is_child_algorithm=True)['OUTPUT']

output_lyr = QgsVectorLayer(final_output, "output_testing", "ogr")

fields = output_lyr.fields()
fieldnames = [field.name() for field in fields]
print(fieldnames)

return {"result": final_output }

The algorithm itself is returning the right output, however I'm not able to access the fields during inside the code, the print is returning me a empty list, however my final_output DO have fields.
Someone, could give me a glimpse of what I've missed?


texts: Converting Excel Data to Projected Point Shapefile says Invalid SQL statement was used?

texts:
I am trying to take an excel file, display xy data, export, and then project. When using the VTDeerCulverts excel file (field names listed in pdf below) the process works until I project the data. At this point the data appears when I identify points but when I try to display the attribute table the following message appears.

After this happened I tried to edit the field names again, as you can see in VDCulvUp (field names listed in pdf below) However when I did this, the problem appeared when I would export the data to a shapefile. At this point the table is empty, no points or data, even though it displays just fine prior when displaying xy data.

PDF of Field Names for Excel Files

Any ideas on a solution?


texts:
Does any one know (OR using OR have experience) exactly what is ArcGIS Marketplace?

According to Google search its related with selling apps and GIS data built using ESRI software's.

Is there any additional information/help?


texts: How to integrate OpenLayers and ArcGIS Server?

texts:
I am using the OpenLayers JavaScript API within a web application.  I am displaying non-base layers hosted on an ArcGIS server.  When I load the ArcGIS layers, I load them using OpenLayers.Layer.ArcGIS93Rest and by specifying the export service URL.  The layer loads perfectly.  However, when the user clicks on a feature, I want to be able to display a popup showing information about the feature that the user clicked on.  Can anyone explain the best way to retrieve feature information?

Thanks in advance!


texts: How to divide a certain country on a map into grids?

texts:
I want to divide a certain country into grids 1,2,3,..., which are divided equally - in other words, have the same width and length. Then, I want to know into which grid a point coordinate will fall and put a mark into that grid. The grids lines will not be visible. For example, if I have points a,b,c, I need to put these points into the grids according to their coordinates. Accuracy is not an issue. I prefer a way that does not require any paid or rate limited services if possible. My application programming languages are Python and JavaScript. I'm using OpenLayers to draw markers on the map.

I'm totally new to GIS systems and maps. I need something to start with. Any suggestion is appreciated.


texts: Using Model to define a group/range of numbers then populate a new field with the new definition

texts:
In one single table I have a rank in one field, and let's say I have rank 1 through 3 in one field and I want to define that group of rankings of 1,2 and 3 as having a score of '10', then ranks 4 through 7 as having a score of'9' and so on, how would I do that in ModelBuilder?

In other words I am trying to define groups of ranks. So far I am trying to look at using python script in a calculate field variable in ModelBuilder.


texts:
I need to select a specific layer in the layer tree before I run my python script.
I have tried
Activating layer by its name in PyQGIS?
 with the code

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

but on the layer tree, there is no change of focus. The layer does not get selected. How can I accomplish that?


texts: extjs 6.2 throws TypeError: a is null and TypeError: contentEl is null

texts:
I am trying to setup my very first working basic example of a GIS using GeoExt 3, Openlayers 4.4.1 and ExtJS 6.2.0 within a Symfony2.8.24 project.

The Symfony part I'm used to work with, I am trying to get the GIS part to work. I have gone through several posts and tutorials just to get this going and when thought I had finally figured out the dependencies I should include in my html file, I'm getting this error I don't even have a clue about. Actually, as you can see in the title, it's two errors, which correspond to ext-all.js and ext-all-debug.js respectively.

In other words:

If I include ext-all.js, then TypeError: a is null shows.

If I include ext-all-debug.js, then TypeError: contentEl is null shows.

Code snippet from my map.js file containing the line which triggers the error:

var description = Ext.create('Ext.panel.Panel', {
            contentEl: 'description',
            title: 'Description',
            height: 200,
            border: false,
            bodyPadding: 5
        });

        Ext.create('Ext.Viewport', {
            layout: 'border',
            items: [
                mapPanel,
                {
                    xtype: 'panel',
                    region: 'west',
                    width: 400,
                    layout: {
                        type: 'vbox',
                        align: 'stretch'
                    },
                    items: [
                        treePanel,
                        description // <-- THIS IS THE TROUBLESOME ONE
                    ]
                }
            ]
        });

And this is the code snippet in ext-all-debug.js where the line with the error resides:

if (me.contentEl) {
    pre = Ext.baseCSSPrefix;
    hidden = pre + 'hidden-';
    contentEl = me.contentEl = Ext.get(me.contentEl);
    contentEl.component = me; // <-- ERROR RESIDES HERE
    contentEl.removeCls([
        pre + 'hidden',
        hidden + 'display',
        hidden + 'offsets'
    ]);
    me.getContentTarget().appendChild(contentEl.dom);
}

Could you tell me why this is happening??


texts: Iterate over feature classes with multiple wildcards

texts:
Is there a way how to iterate over feature classes in Model Builder in ArcGIS 10.3 using a number of different wildcards?

The wildcards I want to use are strings in the middle of the feature class names. If it was just one wildcard, I would use (asterisk)IVORY(asterisk) and it would do the trick but with a number of different ones it doesn’t work. I tried to type “IVORY or CAMEROON” and also “IVORY,CAMEROON” but neither worked.

It would probably be possible to do it creating submodels but when I have a large number of wildcards it would be too much work so I'm looking for a different solution.

Feature class names examples: "feature_class_1_NA300002_IVORY_COAST_AUS15" and "feature_class_2_NA300002_CAMEROON_MAR15" etc.

If possible, I would like to do it without the need to use a scripting language.


texts: Selecting Counties within Region via Spatial Query using QGIS/PostGIS?

texts:
I'm using GQIS 2.2 and PostGIS on Postgresql 9.3.

I have 5 regions in the contiguous United States. Each region has states, and each state has counties.  I would like to think that the data I'm given has each region, state, and county line up cleanly, but one can never be sure.

Here is an example:

I'm at the point in the project where I need to select all counties within a region.  Success is getting all the counties within a given region (say, the Upper Midwest).  None of my attempts at QGIS or PostGIS give me the results I want.  It may either be because my queries are poor (I'm a beginner), or because my data is dirty, or both.

Here's what I've done so far:
I've created PostGIS queries and also tried using QGIS's interface.

ST_Within:

Using QGIS's interface, we can see that ST_Within misses the border counties.

ST_Touches:

Touches seems to fail pretty badly - none of the counties are actually the counties that I want. They seem to be outside the selection! This is what leads me to think that I may have to clean my data somehow.

ST_Intersects:

Intersects is another step in the right direction, but it overselects.

I also tried to get fancy using centroids, but I hit a few interesting edge cases involving centroids falling outside of county polygons and islands.

I could always just use manual selection, but this seems like an interesting problem.  If the problem is just of selection, then fine. But I think the real problem is that my data is off.

I'm seeking any data cleaning techniques. I assume that I'm missing a general technique. I've read How to select by polygon in PostGIS with misaligned layers?, and while the solution was indeed very cunning, it doesn't strike me as solving the underlying problem.


texts:
I am trying to work with a OpenStreet map of Ontario, using the OpenLayers plugin. The problem is that the entire global view renders, slowing down my work.

A couple of solutions come to mind:

Render only the extent of the Ontario shapefile
Limit my ability to zoom out past a view of Ontario, but keep the ability to zoom in. I don't know if this can be done at a project level, or at the layer level.

1 would be the ideal solution. For 2, I've discovered that I can limit the minimum scale visibility, but this brings in other complications:

It lets me zoom out 'past' the visible layer extent. What I would prefer is a hard stop, where I would be unable to scroll out farther.
It lets me scroll up/down/left/right of the view of Ontario. Is this possible to limit?


texts: How can I download OSM maps at fixed scale?

texts:
Is it not possible to download my area of interest (AoI) in my required scale?

I would like to download the AoI in tiles with my specified scale ( SVG or Vector PDF Format).and then  stitch all the tiles in Adobe Illustrator to have complete map of my AoI in one uniform scale.


texts: How to import kml file into R and create a single shapefile?

texts:
I am currently trying to import a batch of kml files in R and to create a single shapefile. Basically, it seems to work fine using st_read and st_write.

But if kml files imported have points lines and polygons together, what's the correct st_write options? I tried to use st_write(a3, "a3.shp",layer_options = c("SHPT=POINT")) (or POLYGON, MULTIPATCH,...), but none of these worked. The error message shows "fail to create feature....".


texts: Interpreting semivariogram with high nugget effect?

texts:
I made a semivariogram in R using package gstat, variogram() function. I want to see if there is spatial autocorrelation in the residuals of my model (species abundance as a function of habitat, across sites spaced a few km to 900 km apart, using a glmm).

My units are in km, and so my interpretation is that the range is just over 100 km until spatial autocorrelation is no longer an "issue". I am wondering if someone can explain why the nugget seems so high? Does this mean that even at similar locations, there is still a relatively high degree of difference? Or, does this wavy variogram mean I should adjust my number of lags and lag distance until I get a more typical shape?

In order to investigate a bit further, I also used the function variog() in package geoR, and used breaks=seq(0,100,10), to try and look just at the closer distances (using same points and same model residuals). This one indicates that the closest points are more different, which also doesn't make sense. Maybe this indicates there is no spatial autocorrelation, and that my model accounts for this already.

I found this excellent source, "Geostats without tears", and on page 51 has some good advice on fitting variograms. By this advice, my first one seems to have the correct range. So this goes back to the first question - how do I interpret this?


texts: Get folder structure in code from arcgis/rest/services

texts:
I have never worked with maps before, and I got assigned on this project with ArcGIS. Noone to ask, cause noone has worked here either so I am turning to you guys. I have to create a tree of checkboxes, and clicking them turns the visibility of the clicked tile on or off. I have managed to find some methods, ArcGISTiledMapServiceLayer and ArcGISDynamicMapServiceLayer. but they don't help me because I cannot give them a link like: http://sampleserver3.arcgisonline.com/ArcGIS/rest/services. but what I need is to take all the folder names from there,

 make them the parents in the tree, then get the services that interest me (in this case there are 3 names in the first folder, don't know which one would be a tile)

 and put them as leaves in the tree.
An example of the behavior I need to reproduce is this:

 in the link I pasted, go to BloomfieldHillsMichigan ,then BloomfieldHillsMichigan/Parcels, click on ArcGIS.com Map from the "View In:" secition, then the "Details" tab, then "Content". And there is the checkbox I need to reproduce, and the text is "Parcels". And all of this has to be in a tree structure, with all the folders we had in the folder section. Is there a way I could accomplish this? Thanks


texts: 1m DEM mosaicking results in a very large file; how can I reduce the size of the output file without compromising output quality?

texts: I am trying to mosaic 792 DEM of 1m resolution using this tutorial (https://automating-gis-processes.github.io/CSC18/lessons/L6/raster-mosaic.html), however it results in a large file (27.9GB) and I get this error "MemoryError: Unable to allocate 27.9 GiB for an array with shape (1, 170000, 44000) and data type float32". I understand that this error is due to overcommit issue but how can I reduce the size of output image initially so that I wont face issue in further analysis like MrVBF etc.


texts:
I would like to install GDAL==2.0.1 on Ubuntu. A lot of possibilities are found on SO, but none did work for me. I'll give an overview of a selection of my attempts:
1. Using pip directly
pip install gdal==2.0.1

returns:

Command "python setup.py egg_info" failed with error code 1 in /tmp/pip-build-xob38erk/gdal/

Note: ez_setup and setup_tools are installed
2. easy_install
easy_install GDAL==2.0.1

returns:

setuptools.sandbox.UnpickleableException:
gdal_config_error(FileNotFoundError(2, "No such file or directory: 'gdal-config'"),)

3. Build from source
Set variables
sudo apt-get install build-essential python-all-dev
export BASE=~
export GDALVERSION='2.0.1'

derive a few others
export GDALFILE="gdal-${GDALVERSION}.tar.gz "
export GDALURL="http://download.osgeo.org/gdal/${GDALVERSION}/${GDALFILE}"
export GDALDIR="$BASE/gdal-${GDALVERSION}"

Download, install and compile for Python
cd $BASE
wget $GDALURL
tar -xvzf $GDALFILE
cd $GDALDIR
./configure --with-python
make
sudo make install

I get the same error as described by Jeremy G.
4. install Python-Gdal
sudo apt-get install gdal-bin python-gdal

This runs with success, but no package gdal is known (with import gdal in python shell or pip freeze).


texts: 'Couldn't load PyQGIS' after upgrade to Raring Ringtail?

texts:
I upgraded my Ubuntu to Raring Ringtail.

After that, when starting QGIS I receive a message:

Couldn't load PyQGIS.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
RuntimeError: the sip module implements API v9.0 to v9.2 but the qgis.core module requires API v8.1

Python version:
2.7.4 (default, Apr 19 2013, 18:30:41)
[GCC 4.7.3]

QGIS version:
1.7.5-Wroclaw 'Wroclaw', exported

Python path:
['/usr/share/qgis/python', '/home/rdk/.qgis/python', '/home/rdk/.qgis/python/plugins', '/usr/share/qgis/python/plugins', '/usr/local/lib/python2.7/dist-packages/pyspatialite-2.6.1-py2.7-linux-x86_64.egg', '/usr/lib/pymodules/python2.7', '/usr/local/lib/python2.7/dist-packages/networkx-1.5-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/psycopg2-2.4.2-py2.7-linux-x86_64.egg', '/usr/local/lib/python2.7/dist-packages/twittercollect-0.1.0-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/pip-1.1-py2.7.egg', '/usr/local/lib/python2.7/dist-packages/Shapely-1.2.14-py2.7.egg', '/usr/lib/python2.7', '/usr/lib/python2.7/plat-x86_64-linux-gnu', '/usr/lib/python2.7/lib-tk', '/usr/lib/python2.7/lib-old', '/usr/lib/python2.7/lib-dynload', '/usr/local/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages', '/usr/lib/python2.7/dist-packages/PILcompat', '/usr/lib/python2.7/dist-packages/gst-0.10', '/usr/lib/python2.7/dist-packages/gtk-2.0', '/usr/lib/python2.7/dist-packages/ubuntu-sso-client', '/usr/lib/python2.7/dist-packages/ubuntuone-client', '/usr/lib/python2.7/dist-packages/ubuntuone-control-panel', '/usr/lib/python2.7/dist-packages/ubuntuone-couch', '/usr/lib/python2.7/dist-packages/ubuntuone-storage-protocol', '/usr/lib/python2.7/dist-packages/wx-2.8-gtk2-unicode']

What is causing such problem? How could I fix it?


texts: Gap filling for linear features using ArcGIS Desktop?

texts:
I have a Landsat time series from which I derived NBR using bands 4 and 7, and delta NBR (to show changes). I then set a threshold based on the mean and standard deviation delta NBR values to reclassify rasters to change (1) and no change (0). The result is a 28 series of images representing change events between years for a study area of about 5,000 km^2. This was all completed in R, and I have the output in raster and vector (polygons) formats. Only pixels labelled as change (1) were used to generate objects (vector polygons).

The output is reasonable for many of the change types occurring in my study area - but I have encountered an issue with roads and other linear features. Instead of the ideal continuous linear (or curvilinear) feature, my output includes many discontinuous features with gaps (as shown in the two examples in the image below).

I have tried to use the Radon transform (PET package in R), but have not been successful. I am looking for a gap filling solution that can be implemented in either ArcGIS or R, using either raster or vector. If the solution is also good for noise removal, all the better! Next steps after gap filling are noise reduction, then classification using Random Forests.

(BTW, this is for an undergraduate project so my level of expertise is moderate)


texts: Plotting points along line at specific distance value using PyQGIS?

texts:
I'm trying to figure out a way to plot points along a polyline at distances that I have in a .csv file (not a point layer with valid lat/long).

There seem to be quite a few ways to create points along a line:

Creating points along line according to attribute data using QGIS?
Creating points in specified distance along line using QGIS？

However, I have not been able to figure out a way based on a given value.

Here is an example of what I am trying to accomplish with this csv data:


texts: Stop QGIS opening in a new workspace in Ubuntu

texts:
Very basic issue I imagine other Ubuntu users will have encountered: every time you open QGIS, it does so in a new workspace - note the arrow in the sidebar is empty not full, just after opening QGIS. I can see how this could be useful, if you want full screen editing and stuff, but it can also be really annoying!

Anyone else had this issues?

Sure there's a very simple solution, I've just not found it yet!


texts: ArcGis Server - GeoData service to serve feature class in json

texts:
I would like to know if there is a way to serve feature class stored in ArcSDE into json format other thant publishing a map service and using queries with the REST API.

For instance, is it possible to publish a geoData service and serve the full features class in json format?

Is there any other way to serve json data from feature class stored in ArcSDE?


texts: LiDAR data in GPS time, only in second. Missing week?

texts: I have LiDAR data LiDAR Data captured, unsure what units or format this time field is in? which I learned that my LiDAR data is in GPS second, however I do not have a column of data which shows the GPS week? Is that something typically collected and should be available in the data?


texts: Load GeoJSON files as a JavaScript variable. Using Require.js

texts:
Trying to load a taluka.geojson

With reference to this post here [this guys post with require][1]

I am sure my GeoJSON file works [Here is a sample][2]

Here is my code:

Here is the error:

require.js:168 Uncaught Error: Module name "taluka.geojson" has not been loaded yet for context: _. Use
require([])
https://requirejs.org/docs/errors.html#notloaded
    at makeError (require.js:168)
    at Object.localRequire [as require]

Here is my file structure just in case:


texts: How to add different sized rasters in GDAL so the result is only in the intersected region

texts:
I'm writing a Python method that adds two rasters and generates a single raster output.  For reasons beyond my control, the extents of the input rasters are different, but they do overlap.

Is it possible to operate exclusively on the input raster pixels that are overlapped in the 2 overlapped regions to generate my ouput such that the output raster extent is only the intersecting region of the two inputs?


texts: Merging two datasets where polygons are nearly identical using geopandas?

texts:
I have two polygon-shapefiles representing similar data. However, the features in these layers are not always 100 % identical. So, a feature in shapefile 1 may be present in shapefile 2, but the form of the polygon is not exactly the same (vertices missing or different). Also, shapefile 1 contains features which are not present in shapefile 2. The opposite is also true.

I'm writing a Python script to merge these datasets. I would like to obtain a layer with all features from shapefile 1, and the attributes of the features of shapefile 2 for those features that are nearly identical to the features of shapefile 1. To define 'Nearly identical' I use the criterium of 75% of the area is intersecting.

How can I efficiently compare the features of shapefile 1 to the features of shapefile 2?

This code works to transfer the attributes of an 'nearly identical' feature of shp2 to shp1. However, it is extremly slow on large datasets, because of the itertuples.

df1 = gpd.read_file('shapefile1.shp')
df2 = gpd.read_file('shapefile2.shp')

df1['attribute1_of_shp2'] = 0
df1['attribute2_of_shp2'] = 0
for feature in df2.itertuples():
    geom = getattr(feature, 'geometry')
    attr1 = getattr(feature, 'attribute1')
    attr2 = getattr(feature, 'attribute2')
    intersection = df1['geometry'].intersection(geom)
    df1['RelIntersection'] = intersection.area/df1.area
    df1.loc[df['RelIntersection'] > .75, 'attribute1_of_shp2'] = attr1
    df1.loc[df['RelIntersection'] > .75, 'attribute2_of_shp2'] = attr2

Is there a way to speed-up this?

I've been looking to geopandas overlay but did not find a working solution yet.


texts:
Having a PostgreSQL database, with the following tables:

CREATE EXTENSION postgis;

DROP TABLE IF EXISTS trips CASCADE;
DROP TABLE IF EXISTS locations CASCADE;

CREATE TABLE IF NOT EXISTS trips (
  trip_id SERIAL PRIMARY KEY,

  start_location bigint NULL,  -- location id
  end_location bigint NULL,   -- location id

  start_date TIMESTAMP WITHOUT TIME ZONE NOT NULL,
  end_date TIMESTAMP WITHOUT TIME ZONE NOT NULL,
  bounds geography(POLYGON, 4326) NOT NULL,
  points geography(LINESTRING, 4326) NOT NULL,
  timestamps TIMESTAMP WITHOUT TIME ZONE[] NOT NULL
);

CREATE TABLE IF NOT EXISTS trips_transportation_modes (
  mode_id SERIAL PRIMARY KEY,
  trip_id SERIAL REFERENCES trips(trip_id) NOT NULL,
  label TEXT NOT NULL,
  start_date TIMESTAMP WITHOUT TIME ZONE NOT NULL,
  end_date TIMESTAMP WITHOUT TIME ZONE NOT NULL,
  start_index INTEGER NOT NULL,
  end_index INTEGER NOT NULL,
  bounds GEOGRAPHY(POLYGON, 4326) NOT NULL
);

CREATE TABLE IF NOT EXISTS locations (
  label TEXT PRIMARY KEY,
  centroid GEOGRAPHY(POINTZ, 4326) NOT NULL,
  point_cluster geography(MULTIPOINT, 4326) NOT NULL
);

CREATE TABLE IF NOT EXISTS stays (
  stay_id SERIAL PRIMARY KEY,
  trip_id SERIAL REFERENCES trips(trip_id),
  location_label TEXT REFERENCES locations(label),
  start_date TIMESTAMP WITHOUT TIME ZONE NOT NULL,
  end_date TIMESTAMP WITHOUT TIME ZONE NOT NULL
);

I am connected to the database and will be making some queries using Flask (Python). This application is supposed to encode the query results in JSON/GeoJSON, to be sent to my front-end.

The problem here is to populate this database. I need to load a folder containing multiple GPX files, parse them and fill these tables with the information I got from the GPX's. Is this attainable with python?

Basically the idea is to get the list of GPX files and parse each one.

Out of the scope of this question, in the future I will convert the parsed results to the formats specified by the tables, and then insert the data in the database using the psycopg library (that I'm already using to connect to the database).

Is there already a good approach to this problem? Or do I need to start it from scratch?


texts: Classical or robust variogram for incorporation into generalized linear model

texts:
I'm modeling counts of organisms over time at eleven locations.  I'd like to account for temporal autocorrelation in the counts, assuming it's present.

As the data are not equally spaced in time, I'm exploring empirical variograms (and comparing variogram model fits) based on the residuals of a generalized linear model assuming independence.

As shown below, the variogram based on the classical semivariance estimator looks quite different from the robust (Cressie) estimator.  Furthermore, a fit of most variogram models to the classical variogram produces singularity errors (with good reason based on the variogram plot) while a decent exponential model (e.g.) can be fitted to the robust variogram.  Unfortunately, it seems as though nlme:::lme (which I'm using for the generalized model) only uses the classical variogram when fitting correlation structures.

Thus, my questions:

(1) Is the use of the robust estimator for constructing the correlation structure justified in this case?
(2) If so, is it appropriate to fit the the robust variogram model in, e.g., gstat or geoR, and then specify/fix the correlation structure in the lme fit?

# Download data
datURL <- "https://www.dropbox.com/s/54q9ocwztt3swap/rsiggeo_variogram.csv"
dat <- repmis::source_data(datURL, sep = ",", header = TRUE)

# Load necessary packages
library(gstat)
library(sp)

# Convert to appropriate class
# Y coordinates are arbitrarily assigned values for each of 11 sites
# (count locations; assumed independent); Y coords were assigned such
# that each site is separated from others by more than the cutoff distance
# used in variogram construction
coordinates(dat) <- ~x+y

# Classic variogram model with example exponential fit
varClassic <- variogram(resid ~ 1, dat, cutoff = 180, width = 180/25)
# Note singularity error; starting values determine 'fit'
expClassic <- fit.variogram(varClassic, model = vgm(1, "Exp", 1, 0))
plot(varClassic, expClassic, xlab = 'Distance (days)')

# Robust (Cressie) variogram model with example exponential fit
varCressie <- variogram(resid ~ 1, dat, cutoff = 180, width = 180/25, cressie=TRUE)
expCressie <- fit.variogram(varCressie, model = vgm(0.3, "Exp", 50, 0.1))
plot(varCressie, expCressie, xlab = 'Distance (days)')


texts: Removing white border of polygons using QGIS?

texts:
I use the graduated color style for the polygons. But there are white borders between polygons. I tried to set Border: transparant, Border stype: No pen, Border width: 0. It did not work.

I found a similar question Removing border of symbol in QGIS? It is about plotting the graduated points.

My current solution is setting the same border color with the polygon manually. But it needs a lot of time. Does anyone know a better solution?


texts: US Geological Survey Standard Projection?

texts: Is there a standard or default projection on maps created by the US Geological Survey (USGS)?


texts: Print Layout Items. setExtent() and attemptResize()

texts:
I'm trying to arrange some map items onto a print layout. Let's start with the map item itself.
Here is my code to add a map item to the print layout. Specifically I have questions about the setExtent() and attemptResize() functions.
In this code below I set the extent with setExtent() before setting the size with attemptResize()
#This creates a new print layout
project = QgsProject.instance()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName('My Layout')
manager.addLayout(layout)

#This is where I configure layout items
map = QgsLayoutItemMap(layout)
rectangle = QgsRectangle(-1350312, -21811, 1741463, 116086) #These are coordinates around my desired map extent
map.setExtent(rectangle)

map.attemptResize(QgsLayoutSize(250, 180, QgsUnitTypes.LayoutMillimeters))

#adds map to print layout
layout.addLayoutItem(map)

This produces the following Print Layout. Notice that the map is "rendering" (Why????). But the sizing of the item is what I want. If I try to export this as a .pdf or image, the final product is blank.

According to this thread, I need to set the map item size before I set the extent. Here is the code, where I have simply switched the order of setExtent() and attemptResize():
#This creates a new print layout
project = QgsProject.instance()
manager = project.layoutManager()
layout = QgsPrintLayout(project)
layout.initializeDefaults()
layout.setName('My Layout11')
manager.addLayout(layout)

map = QgsLayoutItemMap(layout)

map.attemptResize(QgsLayoutSize(250, 180, QgsUnitTypes.LayoutMillimeters))

rectangle = QgsRectangle(-1350312, -21811, 1741463, 116086)
map.setExtent(rectangle)
layout.addLayoutItem(map)

and I get the following print layout. If I zoom in closely, the extent is still correct. But as you can see the image is very small and I have no idea what the attemptResize() function is even doing here, as the units of where the object is on the canvas makes no sense:

What am I not understanding correctly??


texts: I have a lat/long graticule set in my data frame. All looks good, except that the labels on the right and left read horizontally, when I'd like them vertically. I see them on other maps, so I know it can be done, but can't find how to accomplish this in any searches. Ideas?