texts: Create Feature Class From XY Table Using Python

texts:
I have some txt files that are basically attribute tables (object ID, latitude, longitude, area).  In ArcMaps, I have been converting these files into feature classes and saving them in a geodatabase using the following steps:

Right click file

Create Feature Class > From XY Table

Click on the Coordinate System of Input Coordinates button and then
select coordinate system

Save as File and Personal Geodatabase feature classes and output to
selected geodatabase

Since I have a large number of files, I'd like to script this using Python.  I don't have a problem iterating throuygh the file folder and getting a list of all files, but which function from the ARC Python library should I use?  The is Create Feature Class tool the best one I've been able to find so far, but it created an empty feature class.

i have tried opening the Geoprocessing > results tab while creating the feature class manually, so that I can copy as script, but for these actions I am not seeing any entries in the result pane.


texts: What is the difference between InSAR and PolSAR?

texts: How to define those two SAR Technologies?


texts: Generating TIN from features in ArcGIS Desktop?

texts:
I am changing surface in some parts of area I use for project. To do that I draw polygons and add elevation to vertices using "Interpolate shape". Then I create TIN for edited parts of terrain. Next, I change TIN to raster and update terrain with my changes.

Unfortunately the tool "Create TIN" is not generating TIN model for all polygons from the input layer, there are some random gaps in some parts.

Does somebody know what can be the reason of that? Maybe there is a rule I did not hear about.

Edit:

Above you can see polygon which was used for TIN generation, then in the middle is TIN made with the use of polylines and polygons (polygons without capturing Z-value). On the top is another TIN generated just with the use of polygon (Height Field: Shape.Z, SF Type: Soft_Clip). Here looks like the TIN located on the top looks fine, but not all areas look like that.


texts: Testing different values in one CASE WHEN ELSE

texts:
I am making an inventory of certain buildings in a city. In the object attributes, I'm - amongst other infos - linking to pictures of that building, so that in the layout, for each building, you have a map that shows the location of the building, a description of the building and some pictures of it.

There is a possibility to add up to five pictures, but not all buildings will have five (e.g. in some cases, we only have a picture with a general view of the building, in others we will also have pictures with details of certain parts of the construction). First I noticed, that big crosses are shown when a building doesn't have all five pictures, so I added a CASE WHEN ELSE, which makes sure that when there is a picture, the picture is shown, but when there isn't one, there is just nothing instead of a big cross:

CASE
  WHEN "Foto_2" > 0 THEN @project_home  + '/' +   "Foto_2"
  ELSE ' '
END

That works very well, but now, I noticed that portrait pictures are automatically turned and shown as landscape pictures. Is there a way to have an automatic check whether a picture is portrait or landscape and have it adapted accordingly? Because depending on the building, each picture could be landscape or portrait, it's unfortunately not so that for each building picture 1 is landscape, picture 2 is portrait etc.

I would need something like:

CASE
WHEN "Foto_2" < 0 THEN ' '
  ELSE
  WHEN bounds_width($area) > bounds_height($area)
  THEN @project_home  + '/' +   "Foto_2" landscape
  ELSE @project_home  + '/' +   "Foto_2" portrait
END

That doesn't work, but I don't know how to describe the landscape/portrait requirement, nor how to build multiple WHEN ELSE THEN-checks).


texts: Allowing user to download shapefiles via GeoServer

texts:
We provide layers via shapefiles and/or PostGIS PostGreSQL database in geoserver-2.x.x with GCS Geometry.

geoserver-2.x.x itself uses this shapefiles to render maps.

The shapefiles are located on the server machine where geoserver-2.x.x is installed.

I was just thinking to allow user to download specific shapefile(s) of specific layer(s).

I want user to download the shapefile via the website developed using OpenLayers 3 and GeoExt JavaScript libraries with Spring 4,Jsp, Hibernate.

As per this link:

WFS Output for shapefile

Case 1: Layer which is from PostGIS

Case 2: Layer from shapefile (.shp) (located somewhere in D:\shape\Layer.shp on server)

In both above cases will same code work for allowing user to download shapefiles?
If the files may too large or the data table used in PostGIS is too wide, still it will allow user to download the shapefile?
For case 2 will it create a shapefile by itself?
Because my layers are mixing of above two cases. Would I have same logic?
I want to understand how it works. Because I want to give bulk download also for group layers.
Is this possible? I also want to give the same feature in one of my project in world wind.

Any help is invited.


texts: QGIS crashes every time I try to select the labels tab

texts:
Every time I open the Properties box and try to select the "Labels" tab, for any layer in any project, QGIS crashes "unexpectedly" as it says in its Problem Report. I'm using version 2.14.3.

I have tried uninstalling and reinstalling QGIS, but this doesn't help. The labels that I had turned on previously are still displayed, but I cannot edit them. I thought that it could have been a couple of corrupted layer files, but even if I begin a new blank project and create a brand new shapefile, as soon as I try to select the "Labels" tab, the whole program crashes.

Any advice on how to fix this?


texts: Reducing hourly data to daily data and making a timeseries chart

texts:
I'm using GSMaP image collection in Google Earth Engine. Below you can find a script which helps to reduce hourly data to daily data and extract in point scale. By this using this script I would like to make a time series chart, which will demonstrate precipitation events in particular location in given period. I have tried to make a chart, but it did not give me expected graph, as the Chart function is making the according to rows. Below you can also find the received graph. How can I solve this issue?

// Points
var point = /* color: #d63000 */ee.Geometry.Point([69.15, 54.8331]);

// Collection of point
var pts = ee.FeatureCollection(ee.List([ee.Feature(point)]));

// Start and End Dates
var inidate = ee.Date.fromYMD(2018,5,1)
var enddate = ee.Date.fromYMD(2018,5,10)

// Difference between start and end in days
var difdate = enddate.difference(inidate, 'day')

// Import GSMaP data
var gsmap = ee.ImageCollection('JAXA/GPM_L3/GSMaP/v6/operational')
.filterDate(inidate, enddate)
.select('hourlyPrecipRateGC');

// Time lapse
var lapse = ee.List.sequence(0, difdate.subtract(1))
var inidate = ee.Date('2018-05-01')
var listdates = lapse.map(function(day){
  return inidate.advance(day, 'day')
})

// Iterate over the list of dates
var newft = ee.FeatureCollection(listdates.iterate(function(day, ft) {
  // Cast
  ft = ee.FeatureCollection(ft)
  day = ee.Date(day)

  // Filter the collection in one day (24 images)
  var day_collection = gsmap.filterDate(day, day.advance(1, 'day'))

  // Get the sum of all 24 images into one Image
  var sum = ee.Image(day_collection.sum())

  // map over the ft to set a property over each feature
  var temp_ft = ft.map(function(feat){
    var geom = feat.geometry()
    var value = sum.reduceRegion(ee.Reducer.first(), geom, 1000).get('hourlyPrecipRateGC')
    return feat.set(day.format('YYYY-MM-dd'), value)
  })

  // Return the FeatureCollection with the new properties set
  return temp_ft
}, pts))

//Making a chart
var chart = ui.Chart.feature.byFeature(newft);
print(chart);


texts: Creating wedge buffers from many to one in QGIS

texts:
I have two point layers. One point layer has many random points which are start points. The second point layer has only one point which is the endpoint.
I want to create wedge buffers from the start points layer to the endpoint layer.
Something like this:

All buffers have same parameters only angle degree is different


texts: How to display NoData of a raster in QGIS

texts:
How can I show raster pixels with NoData value in QGIS?

In ArcMap there is a button to choose the color for the NoData value. Is there something like this in QGIS?


texts: Properter to get travel time as cost for network analysis

texts:
What I want to do is calculate travel times between points within a road network. I know there are options with GRASS (via v.net) and pgRouting (and possibly many more). I intent to stick to the standard QGIS network analysis though.

I studied all the necessary material in the Network Analysis Library (pyQGIS cookbook) and I am able to obtain the distance in metres for any two points in my network via the standard QgsDistanceArcProperter() and dijkstra's algorithm. What I am not capable of doing is creating a properter that uses a different cost, namely the travel time per segment using the information on speed from its attributes.

Using python QGis network analysis library to calculate short path based on speed and distance? did a good job in translating the speedproperter that is used in the RoadGraph plugin. I am not able to get neither of their properters to run though. QGIS is crashing, so far I couldn't figure out why.

Does anyone know of another already existing properter somewhere that works? Or can point me to something that might bring me one step farther?

EDIT:
It seems that with QGIS3.0 the whole network analysis is revised. Some methods get renamed (e.g. properter -> strategy).
The method is going to be QgsNetworkSpeedStrategy().
The constructor is explained here, can anyone help to translate that into a python class that I could implement within my QGIS 2.18 python script?


texts: Invalid x dimension error in spatial point interpolation?

texts:
I have a point shapefile with altitude values. I have used gdal/ogr to open shapefile and read points(x,y) and altitude(z) values.

It showing error when I am trying to interpolate with scipy griddata modue.

I return back following error:
raise ValueError("invalid number of dimensions in xi")

Please take a look of my codes:

# imports
    from osgeo import gdal, ogr, osr
    import os
    import matplotlib.pyplot as plt
    import scipy.interpolate as il
    import numpy as np
    from itertools import cycle

    # obtaining xyz to interpolate
    a=[]        # to store all x valuesa
    b=[]        # to store all corresponding y values
    c=[]        # to store all corresponding z values
    # input output
    shapefile = r"P:\May2014\Test\LielaisIV_Bp_Schichten.shp"
    outras=r"P:\Interpolatio_spline_with_bariers\outras.tif"
    ## opening a shape file
    DriverName = ogr.GetDriverByName("ESRI Shapefile")
    dataSource = DriverName.Open(shapefile, 0)
    print dataSource
    layer = dataSource.GetLayer()
    srs= layer.GetSpatialRef()
    #looping through features to get geometry reference for each ........
    for i in range(0, layer.GetFeatureCount()):
        # Get the input Feature
        inFeature = layer.GetFeature(i)
        geom = inFeature.GetGeometryRef()
        x= geom.GetX()
        a.append(x)
        y= geom.GetY()
        b.append(y)
        z= inFeature.GetField("of_wth_01")
        c.append(z)
    # obtaining max and min extent
    xmin,xmax,ymin,ymax = [min(a),max(a),min(b),max(b)]

    #size of 1 m grid
    nx = (int(xmax - xmin)+1)
    ny = (int(ymax - ymin)+1)
    # Generate a regular grid to interpolate the data.
    xi = np.linspace(xmin, xmax, nx)#, endpoint=False)
    yi = np.linspace(ymin, ymax, ny)
    xi, yi = np.meshgrid(xi, yi)
    # Interpolate the values of z for all points in the rectangular grid
    zi = il.griddata((x, y), c, (xi, yi),method='linear')

Does anyone have any idea why diention error ocuurs?


texts: How to filter out multipart polygons and then buffer?

texts:
I've written a small script trying to accomplish the task described in the question title. It basically looks like:

in_fc = '...' # Some shapefile goes here

# Get the fields (column names) of the attribute table
desc = arcpy.Describe(in_fc)
shape_field = desc.ShapeFieldName

# Initialize an empty polygon list
singlepolys = []

# Open the attribute table of the shapefile
rows = arcpy.SearchCursor(in_fc)
for row in rows:
    poly = row.getValue(shape_field)

    # Check if single part (not multipart)
    if not poly.isMultipart:
        singlepolys.append(poly)

# Print areas of single polygons - just for debugging
for p in singlepolys:
    print p.area

# Buffer code would go here...

The final printing of the areas shows that only multiples of the last polygon have been copied to the list singlepolys. I really don't understand why.

Any help appreciated!

Edit

A clumsy work-around is possible (I've found a faster one but this is closer to the original thought). Here goes:

# ... Skipping initialization

# Open the attribute table of the shapefile
rows = arcpy.SearchCursor(in_fc)
for row in rows:
    poly = row.getValue(shape_field)

    # Check if single part (not multipart)
    if not poly.isMultipart:
        # Buffer this polygon only and add to list
        geom = arcpy.Geometry()
        result = arcpy.Buffer_analysis(poly, geom, '1000 Meters')
        singlepoly_bufs.append(result[0])

# Now union the buffered features
out_geom = arcpy.Geometry()
union = arcpy.Union_analysis(singlepoly_bufs, out_geom)

# ... Skipping the end


texts:
I have the new Macbook Pro M1 with BigSur. I have installed miniforge via new Homebrew (compatible for M1 chips). Then I created a new conda environment and then installed few libraries, one of them was GeoPandas (conda install geopandas). However, it seemed to failed fetching it and got this error:

UnsatisfiableError: The following specifications were found to be
incompatible with each other:
Output in format: Requested package -> Available versions

Despite the error, I went ahead and created a new Jupyter notebook and wrote only one line "import geopandas as gpd" and did not get any errors. Then I wrote a second line "gpd.read_file('shapefile.shp')" and got this error:
AttributeError: module 'geopandas' has no attribute 'read_file'

In addition, I also tried doing this via pip (before I even installed miniforge but got the same error eventually).
Can I just check, are people using GeoPandas in M1s without issues or am I doing something wrong?


texts: Adding NOAA Radar to JavaScript Application

texts:
I am atempting to add NOAA's Dopler Radar to my web. For whatever reason the later will not display on the map. No errors are being reported and I have used console.log() to confirm the data is being pulled correctly form noaa's arcGIS server. Any thoughts as to how to configure this so the data gets displayed?

maps.addWeather = function addWeatehr () {
    console.log("addWeather Called...");
    var weatherServiceURL = "http://gis.srh.noaa.gov/arcgis/rest/services/RIDGERadar/MapServer"
    var weatherLayer = new FeatureLayer(weatherServiceURL,{
        outFields:["*"],
        visible:true
    });

    console.log(weatherLayer);

    maps.map.addLayer(weatherLayer);
};


texts: How to load layers from postgis to dotspatial or sharemap viewer

texts:
I want to develop a PostGIS viewer for PostgreSQL (with SQL query editor) to visualize PostGIS layers.

Which spatial library (.net C# supported) do I have to use?


texts: How do I update attribute values in a composer attribute table?

texts:
QGIS 2.18.13/Linux Mint 18.1, XFCE (but really Mate) 64 bit.

I have a configured an attribute table in a print composer, it pulls values from an imported (GE) KML layer, "name" and "description" (and lat/lon).  These display nicely.  But I updated the values in GE and these updates are not reflected in the attribute table.  I've pressed the "refresh table data" button and the "refresh view" button (many times) and even restarted to no avail.

Looking at the top values in for the field values in the Print Composer->Attributes->Expression Dialog window,  I get this (old, not refreshed) list:

While looking at the same data via the label tool in map mode Map->Layer->Properties->Label With->Expression Dialog (not composer mode), I see the correct, current values:

I'd think the table data should update both in map mode and in composer mode, but barring that, is there a way to update the attribute table in composer?  Do I need to delete all tables and recreate them every time the underlying field values change to get them to update?


texts: Adding an "m" to the end of every number in a row

texts: What function do i need to use to add an space and an m after every item in a row, using the field calculator in the attributes table. I have a load of numbers (as strings) and would like to add m after each one


texts: WMS GetFeatureInfo not returning all coordinates of geometry

texts:
I have some geoserver layers created from PostGIS tables. I load them on map using OpenLayers and use WMS GetFeatureInfo to get features on mouse click.

I have observed that GetFeatureInfo request optimizes the returned Geometries by merging vertices which are very close. However, for editing, I need to get all the vertices.

Is there a way to tell WMS GetFeatureInfo to retrieve all vertices?

Here is a sample request I am sending:

http://localhost:8080/geoserver/MyWorkSpace/wms
  ?LAYERS=MyLayer
  &QUERY_LAYERS=MyLayer
  &STYLES=MyLayer_Style
  &SERVICE=WMS
  &VERSION=1.1.1
  &REQUEST=GetFeatureInfo
  &BBOX=10064619.024143,2730298.878509,10064682.024975,2730340.381427
  &FEATURE_COUNT=10
  &HEIGHT=556&WIDTH=844
  &INFO_FORMAT=application/vnd.ogc.gml
  &SRS=EPSG:900913
  &X=302
  &Y=265

As commented by user30184, I could use WFS. But it appears that WFS does not take pixel coordinates as input. So I will have to translate the pixel to coordinates, which can be done.
However, in that case how can I consider style of the layer. For example, for a point layer with large radius, user might not be able to select the points even if clicking inside the point. Am I missing something?


texts:
I'm trying to create a local repository to make easier the distribution of plugins within my company.
Although eveything goes smoothly in QGIS 2.14, it seems I can't repeat the same process with QGIS 3.1 !
The structure of my plugins.xml file looks as follows:

First case: downloading a plugin from gitlab

<?xml version="1.0"?>
<plugins>
 <pyqgis_plugin name="plugin_name" version="1.0">
     <description>description here</description>
     <qgis_minimum_version>3.0</qgis_minimum_version>
     <file_name>plugin_name.zip</file_name>
     <author_name>Wanderzen</author_name>
     <download_url>https://gitlab.com/username/plugin_v1/-/blob/master/plugin_name.zip</download_url>
 </pyqgis_plugin>
</plugins>

Second case: downloading a plugin from a network drive

<?xml version="1.0"?>
<plugins>
 <pyqgis_plugin name="plugin_name" version="1.0">
     <description>description here</description>
     <qgis_minimum_version>3.0</qgis_minimum_version>
     <file_name>plugin_name.zip</file_name>
     <author_name>Wanderzen</author_name>
     <download_url>S:\Users\Wanderzen\shared_folder\plugin_name.zip</download_url>
 </pyqgis_plugin>
</plugins>

In both cases the connection to the repository is successful but when I go to the "uninstalled plugins" tab and I click on the plugin I want to install I get those errors:

with Gitlab: Can't install the plugin: failure in the unzipping of the plugin. The plugin may be corrupted or missing in the repository (Échec de l'installation: Échec de la décompression de l'extension. Soit elle est corrompue, soit elle est manquante sur le dépôt)

with the local directory:Can't install the plugin: unknown "s" protocol (Échec de l'installation: Le protocole "s" est inconnu)

Has someone faced the same problem ?!


texts: Cannot access existing worksheet with xlsxwriter using arcpy tool

texts:
I want to populate an existing Excel sheet with values from an ArcGIS feature class. I have made a tool using arcpy and xlsxwriter. I use shutil to copy the template excel in the new location and this is how I create the workbook and the worksheet object:

shutil.copy2(template_xl,out_file)
wb=xlsxwriter.Workbook(out_file)
arcpy.AddMessage(type(wb))
ws=wb.get_worksheet_by_name("POS")
arcpy.AddMessage(type(ws))
ws.write(4,0,"test")

And I get the messages:

<class 'xlsxwriter.workbook.Workbook'>
<type 'NoneType'>
Failed script POSCreator...

Traceback (most recent call last):
  ----, line 62, in <module>
    ws.write(4,0,"test")
AttributeError: 'NoneType' object has no attribute 'write'

If I use:

ws=wb.worksheets()[0]

Then I get the message list out of index, although I need the first (and only) sheet in the excel file.
The variables out_file and template_xl are valid paths.


texts: Need to convert from ISpatialReference to SpatialReference

texts: I am attempting to write a method in C# that will accept as input, an ESRI.ArcGIS.Geometry.ISpatialReference object and convert it to an ESRI.ArcGIS.Client.Geometry.SpatialReference object.  I need this for converting an IPoint to a MapPoint.  Are there any online resources for this?  Thanks!


texts: ArcMap Point Distance giving decimal degrees even when in projected coordinate system

texts:
I'm trying to calculate pairwise river distance among 60 XY coordinates in a small region of Mexico with ArcMap 10.2.2.  I have defined my projections and projected both my stream network and XY coordinate layers as projected coordinates using the NAD_1983_UTM_Zone_14N since the entire area that I am looking at falls within zone 14 as far as I can tell.  I followed all the steps to start a new feature data set, imported my river network and XY layer, built a new network data set, etc.  All of that seems to be fine.  But when I use the Point Distance tool, specifying my XY layer as the input field and the near objects (since I want the distance between all possible combinations of my points), it always gives me decimal degrees.  I also changed the Display setting in the Layer Properties to Kilometers, but that didn't help.  I've tried starting from scratch with a blank map, and projecting everything in projected coordinates, but still no luck.

How do I get Point Distance to give KM or miles rather than decimal degrees?


texts: Resample raster to expand extent in QGIS

texts:
I am building a groundwater model and I need to assign elevations to the base of geologic units (where the unit is present) and at a minimum thickness where the unit is not present. I have a topography raster that covers my entire model area and raster of the geologic units that cover only the area where the unit is present.
I want to use raster calculator to calculate layers base elevation based on whether the unit is present or not.
I tried using a Boolean expression in raster calculator.
For example, for layer 1. I want to assign the base elevation of the layer to the base of FILL where the FILL unit is present and where the unit is not present, assign base of layer to 1 m below the topography.
I was using the following expression
((("topography" - "FILL" >=1)* "FILL") + (("topography"-"FILL"<1)* topography_minus_1)

For some reason, this is not working and it only calculates a really tiny area instead of completely across the domain.
I thought it might have to do with the FILL unit being on a different grid to topography, so I tried r.resample to try and expand the extent of the FILL unit but when I put in the extent of the topography raster and try and resample, it crashes and then QGIS shuts down.
FILL dimensions X: 865 Y: 559 Bands: 1
TOPOGRAPHY dimensions X: 3295 Y: 2056 Bands: 1
How can I make these layers from my rasters based on a boolean expression?


texts: Combining st_snaptogrid with st_centroid and group by using PostGIS?

texts:
I have a table of devices, each of them has a field 'location' that is a postGis geometry(Point,4326)

I need to operate a GROUP BY based of the placement on a grid, I'd use ST_SnapToGrid and then I need to get the ST_Centroid on each set of points in the cell.

I can correctly group by grid:

SELECT count(*) as count,
   st_astext(ST_SnapToGrid(location, 50)) as grid
FROM device_device GROUP BY grid;
 count |       grid
-------+------------------
  3864 | POINT(0 0)
  1462 | POINT(0 -100)
  1377 | POINT(0 100)

and I can get the centroid out of many points:

testgis=# SELECT ST_AsText(ST_Centroid(
            ST_Collect( ARRAY(SELECT location FROM device_device) )));
                   st_astext
-----------------------------------------------
 POINT(-0.281948442456614 0.00443529565569083)
(1 riga)

I cannot understand how to mix the 2 and calculate the centroid for each cell:

testgis=# select count(*) as count,
  st_astext(ST_SnapToGrid(location, 50)) as grid,
  st_astext(ST_Centroid(ST_Collect(ARRAY(location))))
FROM device_device group by grid;
ERRORE:  errore di sintassi a o presso "location"
RIGA 1: ...) as grid, st_astext(ST_Centroid(ST_Collect(ARRAY(location))...


texts: How to extracts points in polygon in GRASS?

texts:
I have a point vector with attributes (cat, dem_value, ndvi_value).

I have a another vector/polygon (converted from contour) with attributes (cat, ID, dem_value). I need to extract points (points with corresponding ndvi_value and dem_value) that are inside the polygon.
I've tried using v.extract, v.select and v.vect.stats. However, I could not get the results.

I used following command in GRASS (wanted to delete the all attributes that are not within the polygon):
v.edit -r map=DEM_exp_point tool=delete polygon=688577.298919,5087607.23962, 688576.065972,5087617.45547, 688582.406842,5087621.33044, 688551.232139,5087597.13671, 688550.880526,5087597.13671, 688562.659564,5087604.69639, 688575.493441,5087606.80606
Selecting features...
619 of 792 features selected from vector map
619 features deleted
Building topology for vector map ...
Registering primitives...
173 primitives registered
173 vertices registered
Building areas...
0 areas built
0 isles built
Attaching islands...
Attaching centroids...
Number of nodes: 0
Number of primitives: 173
Number of points: 173
Number of lines: 0
Number of boundaries: 0
Number of centroids: 0
Number of areas: 0
Number of isles: 0

v.edit complete.

seems it works .. I can see the points (graphical Map display) that lie inside the polygon but no change in the point vector (attribute table). How can I distinguish and extract those points from the Attribute table?


texts: ArcGIS Desktop and Pro read incorectly one band raster as two

texts:
This seems silly, but I have spent too much time searching and making my life over difficult -meaning I am doing something wrong:
I have a TIFF file that is read correctly as one band by tools such as Rasterio, and other Python libraries, but both ArcGIS Desktop and Pro read the file as two-band (R+G), why!? Do I need to explicit request to import the file as one band? And if so how?
Also, is there something missing in the file that makes Esri tools miss the fact that the file is one band?


texts: pgrouting extension slow in PostgreSQL 10?

texts:
I am trying to implement routing on my road dataset using pgrouting function pgr_createTopology. It works fine on postgresql 9.6  but gets extremely slow and get stuck as nodes increase in postgresql 10. Both postgresql 9.6 and postgresql 10 are on the same computer.

This issue occurs when the dataset is large. For small datasets this function works completely fine in both postgresql 9.6 and 10.

I can not understand the reason why pgr_createTopology functions works completely fine in postgresql 9.6 and gets stuck in postgresql 10.

Is there anyone who could guide me regarding this issue?

(postgis version is 2.4.4 and pg_routing version is 2.6.0)


texts: Creating a custom CRS for the South African coordinate system (SACRS) in QGIS?

texts:
My question relates to creating a custom CRS of the official coordinate system used in South Africa, called SACRS.

The SACRS is a projected coordinate system that must be referenced to BOTH the "Hartebeesthoek94" datum and the "Gauss Conform Coordinate System" (see help guides at the bottom of this post). I seek technical advice on creating this projected coordinate system as a custom CRS.

What should I type in the parameters box in QGIS (proj4 format)?

I am not sure if any of the following coordinate systems are already close (or similar) to the SACRS, but perhaps they are a good starting point for a custom CRS:

Hartebeesthoek94 / Lo25 [EPSG 2051]
South African CRS: HBK_NO_25 [ZANGI:ZANGI:HBKNO25]

Background information on the SACRS is available at:

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=13&cad=rja&uact=8&ved=0ahUKEwizmdnHwoXSAhVpLMAKHY3CATUQFghRMAw&url=http%3A%2F%2Fwww.ngi.gov.za%2Findex.php%2Fresources%3Fdownload%3D37%3Athe-south-african-coordinate-reference-system&usg=AFQjCNFywI3HdttJQWFzPXFqtAr4leRdiQ&sig2=mGJdK_hyvRsJXU8-F6w86g

https://ace-notebook.com/An-introduction-to-coordinate-systems-in-south-africa-free-related-pdf.html


texts: DroneDeploy generated Ortho shows skewed in QGIS (projections OK?)

texts:
I've recently bought a drone and have been having a lot of fun enhancing maps with very detailed orthomosaics.

The images place correctly once loaded in QGis, absolute positioning is good enough for what I want to do etc.

The images generated should be WGS84, in Qgis, when looking at the properties that is also reflected. I've set the project properties based on the image.

I've got a very nice and square tennis court, when I view it on the DroneDeploy website it shows as a proper square. However, once loaded in QGIS it's pretty badly skewed.

I'm a (Q)GIS newby and rely on Google to learn things.

I've not been able to figure out what is going on past the checking of projections.

Any suggestions, it's a private property so I can't link to the full image. But if an image is useful I can add some partial screenshots


texts: Playing with Events of BaseTool to draw Geometries.in ArcGis using ArcObjects C#

texts:
I have been trying to create a tool that creates a line using points. The idea in my mind is to get the points from a onmousedown Event and putting them in generic list and draw the line  geometries using those points onmousemove exactly like placing a  transmedia  from a Model Pallate in ArcMap.

I am able to place points on map. But I am unable to create lines on move. I was thinking if I would get the point on onmousemove but was clueless how to do that. Even how to stop the onmousemove. I do not try to use the line feedback or any handler for creating this. I was thinking to create the tool from scratch not taking help from any other tool.


texts: Splitting irregular polygon into equal areas using QGIS?

texts:
I am a surveyor and use QGIS for maintaining soil test records on civil projects and am constantly provided with test lot area or number of lots required. These are more often than not irregular in shape and am wondering.

Is there a plugin/solution to cut a layer that contains irregular polygons into equal areas or number of lots with a specified area, with the result being the most number of lots based on the area provided?


texts: geopandas dissolve versus postGIS postgresql aggregated groupby dissolve

texts:
I have a geodataframe for which I want to dissolve the geometry based on an attribute (fao_id). In geopandas this is easy as gdf.dissolve(by="fao_id", aggfunc="first") but unfortunately this fails due to some invalid geometries.

I ingested the geodataframe in PostGIS to make use of the st_make_valid functionality that geopandas does not have. This method works like a charm but now I want to dissolve my table using postGIS similar to the gdf.dissolve function.

I found this article but the result is a less parsimonious SQL statement:

CREATE TABLE newtable AS
SELECT MIN(fao_id) as fao_id,
       MIN(attribute1) as attribute1,
       MIN(attribute2) as attribute2,
       ST_Multi(ST_Union(t.geom)) as geom
FROM inputtable As t GROUP BY fao_id

some comments: 1) I have to manually type all my columns in the statement which can be tedious with larger tables. 2) there is no native FIRST() aggregate function but using MIN() our MAX() works since all values per fao_id are identical.

I was wondering if there was a better way. My attributes2 3 etc are only dependent on the composite key fao_id. see below:


texts:
Middle of a project, wondering if there is a better historic topo source than TheNationalMap for finding historical maps for an area.
Coring some lakes, would be nice to know that they've been there for at least a few hundred years.


texts: How would I compute the minimum union bounding box given a list of bounding boxes?

texts:
Constraints:

Assume -90 < latitude <= 90 and -180 < longitude <= 180.

A bounding box has the form [minLon, minLat, maxLon, maxLat].

Concerns:

A solution should be clever to realize that 2 boxes on either side of a boundary should be merged over that boundary.

Ex: Box_A has left/right longitude bounds of (177, 179), and box_B has (-179, -177). The result should compute (177, -177) and NOT (-177, 177).... may have confused this here, but look at the pictures below to see what I am getting at.

If 2 boxes are disjoint, (in my head) there are 4 possible union boxes:

Box A is merged with box B on either the left, right, top, or bottom of box A.

More than 2 boxes could really slow down a naive approach.

Question:

Is there a python library that will facilitate solving the following problem:

Given a set of bounding boxes {B}, compute a minimum bounding box C such that for all b in {B}, b is contained in C

, and if not, is there at least a library that can compute this result given lat/lon pairs manually extracted from corner points of these boxes?

Clarity:

I am not concerned about what metrics are used to solve the "minimum" box, whether it be area, boundary differences, etc. As long as it is a reasonable and efficient solution.

Visual:

Input:

Output (red is good, black and purple are bad):

EDIT: I just realized the purple box does not actually make sense here on the flat map (read Logan's comment below). I'm going to leave it for now since the principle still applies to the black box at least.


texts: Change the processing extent of raster in ArcGIS Pro

texts:
I have two rasters. One shows slope:

and the other shows Euclidean distance:

As you can see, the Euclidean layer is smaller than the slope one.

How do I change the extent of the small layer to the extent of the big layer?

I know than in ArcGIS Desktop I could do it in the properties of the layer and then choose processing extent and same as layer X, but I can't find how to do it in ArcGIS Pro.

My end goal is to have the small distance raster  have the extent of the big raster.


texts: Clipping out polygon into buffer layer so that buffer layer shows dominant feature from polygon layer?

texts:
What I am trying to do is clip out a land use layer for the city of Boulder, CO within the confines of a 500ft buffer that I have created for a bunch of road segments. Beyond that I'm trying to get each road segment buffer to only show the dominant land use type.

Here's an image to help visualize a bit:

I have the road segment buffer selected and you can see the land use polygon underneath. Obviously here the yellow land use type is dominant and, as such, I would like for the buffer to clip that out and apply it to the buffer entirely. Other road segments aren't as clear though so I need some way to compute and automate this.

Is this possible?

If not, I'm just trying to show the dominant land use type is per road segment in order to see how a road, and its users, are impacted by said land use as part of a larger study, but I'm at a loss for how to actually do this. This is what I thought would be best, but I'm open to another method entirely if you have an idea.


texts: Ordering tool display within a Python Toolbox

texts: I have created a python toolbox containing a number of classes which define the various tools within. As these classes are part of a workflow in which the tools are to be executed sequentially using outputs from the previous tools, it would be very handy if the tools were presented in this order within the box, reinforcing the workflow. However, they appear to be displayed alphabetically. Is there a way to define the order? The tools are entered into the class Toolbox attribute list, self.tools[tool1, tool2] in the proper order, but they are subsequently alphabetized before being displayed.


texts: Creating a table from the result of ST_DumpValues()

texts:
I have, I believe, successfully loaded a raster into my PG DB. I see the table in PGAdmin but am unable to see it in QGIS. The table in PG has 1 row. Running, e.g. st_count(rast) on the table returns over 600,000. I am,

Unable to view the dataset in QGIS
Unable to run PG raster function (e.g. ST_DumpValues, ST_Aspect).

There is obviously a first step to do after loading data to the database that I am missing. I cannot find info related to this anywhere. What is it that I am missing?

CREATE TABLE raster1 AS SELECT
 ST_DumpValues(ST_Aspect(rast, 1, '32BF'))
 FROM raster

Gives me,

  ERROR: column "st_dumpvalues" has pseudo-type record


texts: Find the intersection angle between two lines counter-clockwise around a circle ArcGIS 10.x

texts:
I am trying to create a tool for measuring the angle between two intersecting lines. I have created a buffer and clipped the lines to within this buffer and computed the angle between the lines using the formulas below:

dx = !X1a!-!X0a!

dy = !Y1a!-!Y0a!

dxa= !X1b!-!X0b!

dya= !Y1b!-!Y0b!

r  = math.sqrt(math.pow(!dx!, 2) + math.pow(!dy!, 2))

ra = math.sqrt(math.pow(!dxa!, 2) + math.pow(!dya!, 2))

c = math.asin(abs((dxdya - dydxa))/(r*ra)) / math.pi * 180

((source How to find the angle between intersecting features in two featureclasses?, user: whuber))

This returns the angle between the two intersecting lines, but, I need the angle between the lines to always be between line1 counter-clockwise to line2. The tool will be used to calculate the angle between a road intersecting with a rail road.

Does anyone know the formula to do this?

The image below shows a sample of the data I am working with and the angle I would like to calculate. The diameter of the circle (buffer polygon) is 10 centimetres.


texts: how can I change the border of my layer?

texts:
I have 2 layers:

one with neighborhood borders
another with service areas (shapefile)

I would like to make one layer with all the data from ServiceAreaMETRO and for the zone where i have no information in ServiceAreaMETRO but is inside of the layer 200001445 I would like to give value 1.
How can i do it?


texts: Get Raster Values from a Polygon Overlay in Opensource GIS Solutions

texts:
I have two layers. A polygon-shape-layer with many tiles and a raster-layer containing CORINE 2006 land cover with many categories in a colourmap.
I want to obtain for each polygon in the shapelayer a sum of each land cover category of the raster-layer.

For example there is a polygon with id '2' and i want to Attributes like this for this polygon (in percent or square meters):

Arable land : 15 %
Forest: 11 %
Streets:2 %
(... and so one)

I tried to do it in grass, qgis (no function), saga (just sums up every to a total value) r(total sum), but i still found no solution. Most plugins (zonal statistics in qgis) only support 0-1 raster layers. v.rast.stats didn't help either. Iam open to any good and smart solution!. Maybe i even used a wrong approach or made mistakes.

In Arcgis this task is quite easy, if am remember right, but i am still missing a good solution for your everyday linux user.

I am running a debian linux system and this why i can only use programs for this OS.

EDIT:
Because this question still has so many views and visitors:
I wrote a QGIS-plugin, which also is capable of calculating the landcover of raster layer. I have'nt coded a polygon overlay yet, but it definitely planed. Find the plugin here and install the Scipy library first.


texts: How to expedite process of creating features

texts: I am creating TIN out of each two cross sections and in the end I merge them all together and create a seamless TIN using python scripting. First I create all the TIN and save them somewhere in my local hard drive and finally I merge them. Someone told me there should be a way to do this process in the memory to save time and them create the final TIN and copy that in my hard drive. I already use feature layer for cross sections to expedite the process. I was wondering if there is still any trick I can use to lower the process time?


texts: Call Sql Server stored procedure from arcpy

texts:
Is there a way to execute a SQL Server stored procedure with arcpy or Python and retrieve the results returned by the stored procedure?

There is a very similar question posted on ArcGIS forum with no answer so far.


texts:
I was wondering if anyone has used the i.albedo processing tool in QGIS? I am looking to do an analysis on a composite image from Landsat 8 and don't really know where to begin. I am using version 3.12 and keep getting an error when I open up the i.albedo algorithm.
Thanks,
Nick


texts: Is there a way to decrease the auto-hide delay for dockable windows in ArcMap 10?

texts:
I'd like for my dockable windows to hide a little faster than they do now.  Is this a user or application setting within ArcMap or ArcGIS?

More details: I'm specifically interested in ArcMap, but if it's one setting for the entire ArcGIS desktop--I'll take it.


texts: I want to perform atmospheric correction on Landsat 8 images in ENVI 5.1. In order to create radiance image, I selected basic tools> preprocessing> calibration utilities> Landsat calibration. However, OLI sensor does not exist in sensor type. Would you please help me how I can perform this analysis in ENVI 5?


texts: Combo box in python add-in for Arc - keep selected value visible

texts:
I have a python add-in for Arc using the addin-assistant. It populates values in one combo box based on the selection from the other.  All of that works fine.

What I can't figure out is how to keep the selected string from the first box visible. See the image... How do I get the "Zoom to:" combo box to display the value that the user selects?

I must be missing something but I cannot figure out what I am missing or where to put it. Code below.

import arcpy
import pythonaddins

class ComboBoxClass7(object):
    """Implementation for arcpythonAddins_addin.combobox_1 (ComboBox)"""
    def __init__(self):
        self.items = ["quad", "county"]
        self.editable = False
        self.enabled = True
        self.dropdownWidth = 'WWWWWWWW'
        self.width = 'WWWWWWWW'
    def onSelChange(self, selection):
        combobox.enabled = True
        global fcName, fcField
        if selection == "quad":
            fcName = "quads"
            fcField = "quad"
        else:
            fcName = "counties"
            fcField = "COUNTYNM"
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        pass
    def onEnter(self):
        pass
    def refresh(self):
        pass

class ComboBoxClass1(object):
    """Implementation for arcpythonAddins_addin.combobox (ComboBox)"""
    def __init__(self):
        self.editable = True
        self.enabled = False
        self.dropdownWidth = 'WWWWWWWWWWWW'
        self.width = 'WWWWWWWWWWW'
    def onSelChange(self, selection):
        arcpy.SelectLayerByAttribute_management(layer, "NEW_SELECTION",  fcField + "='" + selection + "'")
        self.df = arcpy.mapping.ListDataFrames(self.mxd)[0]
        self.df.zoomToSelectedFeatures()
        arcpy.SelectLayerByAttribute_management(layer, "CLEAR_SELECTION")
        arcpy.RefreshActiveView()
    def onEditChange(self, text):
        pass
    def onFocus(self, focused):
        self.mxd = arcpy.mapping.MapDocument('current')
        global layer
        layer = arcpy.mapping.ListLayers(self.mxd, fcName)[0]
        items = []
        with arcpy.da.SearchCursor(layer, [fcField]) as cursor:
            for row in cursor:
                items.append(row[0])
        self.items = sorted(items)
    def onEnter(self):
        pass
    def refresh(self):
        pass


texts: How did OpenLayers 3 managed to export map from other domain?

texts:
As far as I know, you cannot use canvas.toDataURL() on a tainted canvas (i.e. cross-origin). Looking at this example it seems the tiles are coming from different origin (tile.openstreetmap.org). How did OL3 managed to use toDataURL() in this case?

I tried monitoring network traffic and it seems no proxy are being used to circumvent cross-origin rules.


texts: Setting field defaults/domains based on subtypes

texts:
I have a rather complicated FeatureClass, in which one of the fields has a different default value, and a different domain, depending on the value in the nominated SubType field.

Using the rather useful code in this example, I get this output: (trimmed to just the field in question, for length)

Code: 1
Default: True
Name: EducationalFacility
SubtypeField: FCSUBTYPE
Fields:
 --Field name: FACILITYTYPE
 --Field default value: 3
 --Domain name: dEducationalFacilityType
Code: 2
Default: False
Name: HospitalFacility
SubtypeField: FCSUBTYPE
Fields:
 --Field name: FACILITYTYPE
 --Field default value: 6
 --Domain name: dHospitalFacilityType
Code: 3
Default: False
Name: CommunityFacility
SubtypeField: FCSUBTYPE
Fields:
 --Field name: FACILITYTYPE
 --Field default value: 3
 --Domain name: dCommunityFacilityType

Note the different domains for each subtype value. The FACILITYTYPE field itself has no domain set.

When I copy this FeatureClass, using FeatureClassToFeatureClass with arcpy.env.transferDomains = True, I lose all this interesting subtype information. I can recreate the domains for fields that have domains set directly, but I can't work out how to transfer this domain-within-a-subtype stuff with arcpy.


texts: Can I get Geoserver to return a WMS in gnomonic projection?

texts:
I need to be able to pass in the center lat and long.  I have tried using the AUTO Namespace for SRS (which works) but this only supports a few projections, the closest being orthographic.

I have found this link for ArcGIS http://resources.arcgis.com/en/help/main/10.1/index.html#//015400000605000000 which says I could pass in the whole projection, rather than just an EPSG code.  I haven't actually been able to find if this even possible with GeoServer.

I am using GeoServer 2.5.1.


texts: How convert a genuine UK shapefile(other than the one from natural earth) to d3/topojson?

texts:
I am following Mike's page (http://bost.ocks.org/mike/map/) to make UK's map. I've converted Mr. Crone's UK postcode shapefile into topojson but failed to display it correctly.

http://www.sharegeo.ac.uk/handle/10672/51?show=full

ogr2ogr -f GeoJSON postcodes.json gb_pca_2010.shp
topojson -o uk.json --id-property Label postcodes.json

The resulting json is of 2.2M big. I have used the following code to display it

var height, svg, width;
width = 960;
height = 1160;

svg = d3.select("body").append("svg").attr("width", width).attr("height", height);

d3.json("uk.json", function(error, uk) {
  var projection, subunits;
  if (error) {
    return console.error(error);
  }
  subunits = topojson.feature(uk, uk.objects.postcodes);
  projection = d3.geo.mercator().translate([width / 2, height / 2]);
  return svg.append("path").datum(subunits).attr("d", d3.geo.path().projection(projection));
});

The result is that firefox freezes and in the end, I get half black and half white display. :( I have also tried mapshaper to simplify it to 1MB but still not working.

What else I could have tried?


texts: Extract Labels (city names, state names) with shapefiles using QGIS

texts:
I want to extract shapefiles with labels on it using QGIS but labels cannot be extracted along with shapefiles.

Here's the image:

I want this as my output after extraction:

How to do something like this in QGIS? I want to upload this shapefiles in Google Fusion Table as KML but the Labels cannot be extracted?

Or, can you recommend me another software to do this?


texts: Counting intersections between line and points using Field Calculator?

texts:
I'm using QGIS 3.12 and let's say I have a Layer that consists in Lines that represent 'Rivers', and I have another Layer that consists in Points that represent 'Bridges'.
Can I count the number of intersections between these two Layers using Field Calculator?

In the example given the objective would be to have the attribute 'Inters_Bri' be equal to 5 using the Field Calculator because I have 5 bridges intersecting the Blue River.


texts: QGIS 3.x closestSegmentWithContext returns sqrt distance 0.0 and point doesn't snap to polyline

texts:
I have a point and line shapefile.

I am new to QGIS and Python, but have been able to get a working method to snap some of the points to the line. The method I am using is "closestSegmentWithContext", which works where the points are further from the line, but when the points are close (i.e. less than 5-10 meter) it returns a sqrt distance(?) of 0.0, and cannot find the closest point on the line. See for example the below outputs:

Var1:  (4.4249768536154936e-08, <QgsPointXY: POINT(27.9118715303157181 -26.09893318956183705)>, 35, 1)
Var2:  (0.0, <QgsPointXY: POINT(27.94594999999999985 -26.09562199999999876)>, 65, -1)

With Var1 I am able to snap the point onto the line, but the Var2, the point cannot be snapped onto the line and returns the coordinate of the point.

I am using EPSG:4326 - WSG 84, and I am assuming that the projection is in degrees, although I have set the distance to be in meters in the project settings.

Is there anything that I need to add to my code to get the points to snap onto the line.


texts: What data is sent from GPS satellites to receivers?

texts:
I know more or less how the GPS system works. But I would like to understand it in more details. Exactly what data is sent from a GPS satellite to a GPS receiver?

I think there is multiple types of data, since when I start up a GPS receiver, it takes some time before it becomes useful. As what I know this is because the receiver waits on some data that aren't received very often. I guess that the GPS signal that is used for the position is sent more often, but I don't know.

Is there anyone that can explain this in more detail? Or have a useful link?


texts: Seeking alternative to ArcGIS Water Utility Tools to create lateral lines?

texts:
I need to automatically create perpendicular lines from a point to a line (curb box to water main) and also create perpendicular lines from a point to a line or polygon (curb box to building footprint) in order to create water laterals (line).  I don’t need to model, etc., but some way to create this data without manually drawing in all the lines.

I have ArcGIS-Basic 10.6 (license level not compatible with free water utility network tools), X-Tools, AutoCAD Map 2016 and QGIS software.


texts: Alignment problem in 2 vector shapefiles with the same projection

texts:
I have 2 vector shapefiles of Manila that are supposed to align, but they do not (screenshot attached). I have projected them to the same projection (World Mercator WGS 84, EPSG 3395). One shapefile appears to be shifted to the right and does not align with the Google Street basemap and the other shapefile. I am using QGIS 2.14.16.

The pixelised shapefile in the picture was actually from a MODIS 500m raster image which I polygonised. The other shapefile is the administrative outline of Manila, from GADM.

I understand that MODIS data may have some misalignment issues that may be the cause of this problem, which I have not seen in other forums.

Any ideas as to what could be the problem and solution?

Tried to change the shapefiles to EPSG 3857 to match Google Street's projection as suggested by AndreJ. Unfortunately both shapefiles are more misaligned.


texts: Geotools RPFTOCReader NegativeArraySizeException in reading raster

texts:
System: Eclipse IDE in Windows 10 VM Java 8

My code segment that the exception is being thrown:

private void displayLayers(File rasterFile) throws Exception {
    Hints hints = new Hints();

    reader = new RPFTOCReader(rasterFile, hints);

    final MapContent map = new MapContent();
    map.setTitle("Jonathan's attempt");

    Style rasterStyle = createGreyscaleStyle();
    Layer rasterLayer = new GridReaderLayer(reader, rasterStyle);

    map.addLayer(rasterLayer);

    frame = new JMapFrame(map);
    frame.setSize(800, 600);
    frame.enableStatusBar(true);
    frame.enableToolBar(true);

    JMenuBar menuBar = new JMenuBar();
    frame.setJMenuBar(menuBar);
    menuBar.add(new JMenu("Raster"));

    frame.setVisible(true);
}

the createGreyscaleStyle() method is the same one from this tutorial.

Exception:

    Jan 02, 2018 2:37:41 PM it.geosolutions.imageio.gdalframework.GDALUtilities loadGDAL
    INFO: GDAL Native Library loaded (version: 2.2.1)
    Exception in thread "main" java.lang.NegativeArraySizeException
        at it.geosolutions.imageio.gdalframework.GDALImageReader.readDatasetRaster(GDALImageReader.java:331)
        at it.geosolutions.imageio.gdalframework.GDALImageReader.read(GDALImageReader.java:952)
        at org.geotools.coverageio.RasterLayerResponse.readRaster(RasterLayerResponse.java:440)
        at org.geotools.coverageio.RasterLayerResponse.createCoverage(RasterLayerResponse.java:210)
        at org.geotools.coverageio.RasterLayerResponse.compute(RasterLayerResponse.java:180)
        at org.geotools.coverageio.BaseGridCoverage2DReader.requestCoverage(BaseGridCoverage2DReader.java:606)
        at org.geotools.coverageio.BaseGridCoverage2DReader.read(BaseGridCoverage2DReader.java:387)
        at org.geotools.raster.ImageLab.createGreyscaleStyle(ImageLab.java:189)
        at org.geotools.raster.ImageLab.displayLayers(ImageLab.java:107)
        at org.geotools.raster.ImageLab.getLayersAndDisplay(ImageLab.java:95)
        at org.geotools.raster.ImageLab.main(ImageLab.java:69)

    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr
    OpenGL Warning: crPixelCopy3D:  simply crMemcpy'ing from srcPtr to dstPtr

Tracking the stack trace

The exception is thrown within the createGreyscaleStyle() method within my code (link to method is mentioned above) which is called from the Style object I attempt to make within my code excerpt. Within the createGreyscaleStyle(), the reader.read() function is called which eventually calls readDatasetRaster() (method in GDALImageReader).

On line 327 of GDALImageReader the calculated bufferSize is -1173618688, dstWidth is 75264, dstHeight is 41472, and when pixels is first calculated it's -1173618688.

When I run gdalinfo on the A.TOC file I get a size that is 75264, 41472 so it looks like the reader is extracting accurate data.

Problem:
Don't know how to account for the negative value of bufferSize.


texts: Difference between `system:index` and `system:time_start` properties in Google Earth Engine

texts:
If I want to know the date of image acquisition, would it be the system:index or system:time_start property?

I'm afraid I have been created a mess and used wrong dates.

For example, if I have this list, is the date that is referenced by the red arrow the acquisition date or the publication date?


texts: Migrating geodatabase data into PostGIS without Esri apps?

texts:
Can you tell me how I can load geodatabase data into PostGIS without having to install Esri apps?

The solution I found online talked about using ArcGIS to convert GDB to SHP and then importing SHP to PostGIS.


texts: arcpy CompositeBands for each folder, name new raster by folder

texts:
I am processing multibands remote sensing rasters (i.e. Sentinel data). Each set of bands is stored in individual folder. I would like to list folders and list rasters (bands) in each folder to create a raster composite. Each new raster I would like to store under the folder name. How can I name the new raster by folder name?

# create raster composite per folder

import arcpy
import os

#
arcpy.env.workspace = "C:/Users/clipy/"
outws = "C:/Users/output"

# list folders
folders = arcpy.ListWorkspaces()

arcpy.env.overwriteOutput = True

for folder in folders:
    arcpy.env.workspace = folder
    rasters = arcpy.ListRasters("*.tif")
    name = os.path.join(outws, folder[1])
    arcpy.CompositeBands_management(rasters, name)

print "Raster composite complete"

I'm using ArcGIS 10.4 and python 2.7.


texts: Raster to ascii tool converting all data to one value (ArcMap)

texts:
I am having an issue when converting a raster to ascii. With extract by mask I have selected the area I am interested in and subsequently resampled the data from 5km to 2 by 2 m (required for use in modelling system). I need the raster to be in ascii format and using the raster to ascii seems to change all the data to one value. I have used this methodology numbers of times before and this has never caused this issue before.

I wonder if anyone may have any ideas?


texts: What is the best technique for interpolating soil depth from sample points

texts:
This is my first post here.

I have some points from a soil depth survey and I need to interpolate the depth for a whole area.  What is the most appropriate technique for this and is there anything I need to take into account when using the technique?  A pointer to a tutorial would be brilliant!

Thanks


texts: Right coordinates for data defined position of diagrams

texts:
I'm using QGIS 2.8.2 and want to define the position of diagrams with x and y coordinates.

The same thing works with the labels, but when I define coordinates for the diagrams they just disappear.
I've tried different ways of coordinates but everytime they disappear.
I think there is a problem with the Coordinate Reference System. I had to use WSG 84 EPSG:3857 and EPSG:4326 because otherwise my map is stretched.

I'm searching the web and the documentation for hours but i can't solve the problem, i am still a very beginner. Has anyone an idea to fix my problem?


texts: Multiprocessing ERROR 000210 Cannot create output 'name' Failed to execute Intersect

texts:
I am trying to use multiprocessing to process intersects between polygons and points for a list of states and territories (the data sets for the points are large, so my hope was that the multiprocessing would speed the process up). But there seems to be an issue. Suppose I was just doing DC and PR. The intersect first produce an empty table output for PR and then throws the following error for PR:

Traceback (most recent call last):
  File "C:\Projects\hh_intersects\just_hh_intersects.py", line 68, in <module>
    pool.map(multi_proc, state_list)
    ....
arcgisscripting.ExecuteError: ERROR 000210: Cannot create output
T:\...\civis_hh_0520.gdb\PR_hh_intersects
Failed to execute (Intersect).

I'm not sure what I am doing incorrectly here. Any ideas on how to fix this?

import arcpy
import os
import time
import multiprocessing as mp
from multiprocessing import Pool, cpu_count

# Paths for hh_pts_gdb, input_polys_dir, and output_gdb

state_list = ['DC', 'PR']
start = time.time()
p_count = cpu_count()

# Functions
def multi_proc(state):
    state_start = time.time()
    arcpy.env.workspace = output_gdb

    # Make output based on state
    output = os.path.join(output_gdb, "{}_intersect".format(state))

    # Pts
    pt_name = "{}_hh_pts".format(state)
    hh_pts = os.path.join(hh_pts_gdb, pt_name)

    # Polys
    st_poly = os.path.join(input_polys_dir, "{}/{}_UnionOutput.gdb/{}_Union".format(state, state, state))

    # Intersect selected pts with polys
    arcpy.Intersect_analysis([hh_pts, st_poly], output, "ALL")
    print("Output row count for {}".format(state))
    counter = arcpy.GetCount_management(output)
    print(counter)
    state_end = time.time()
    print("{} Process duration: ".format(state))
    print(state_end - state_start)

if __name__ == '__main__':
    pool = mp.Pool(p_count)
    pool.map(multi_proc, state_list)
    end = time.time()
    print("Total Process duration: ")
    print(end - start)

*Note- I have run the following by itself (outside of the multiprocessing) and it worked with state = "DC": arcpy.Intersect_analysis([hh_pts, st_poly], output, "ALL").


texts: I have download a tiff file with elevation from http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp and I try to export in file text or Excel the coordinates and the elevation and I don't know how to export this data in text file.


texts: Shortest and longest amount of numbers that latitude and longitudinal coordinates may have?

texts:
What is the shortest and longest amount of numbers that latitude and longitudinal coordinates may have?

For instance i know that London is 51.5072 0.1275. So the first value is 6 numbers and the second is 4 (or 5 depending on how you look at it).


texts:
I am new to ArcGIS. I have the EDN license (Desktop and Server); have setup ArcGIS for Server in our offline environment; have published several base maps (with caches that ESRI provided me on a DVD) as services; and can use those services with an ArcGisTiledMapServiceLayer from a Silverlight client we are building.

The map displays, I can add our graphics layers to it, and everything works great. But, I do not want to use the default projection that these base map caches have (the 102100 WKID as seen in this arcgisonline service http://services.arcgisonline.com/arcgis/rest/services/World_Street_Map/MapServer). I'd like to use the projection from here: http://services.arcgisonline.com/arcgis/rest/services/ESRI_StreetMap_World_2D/MapServer which is the WKID 4326.

Given just the base map caches that I have, is there any way for me to reproject the map? Can the client do it on-the-fly? Or can I configure ArcGis for Server to do that reprojection dynamically? Do I need to rebuild the whole cache some how to the new projection? Or find a new cache entirely?


texts: Using SearchCursor and UpdateCursor to get multivalues from one table/field to single value of another table/field?

texts:
I'm trying to use a combination of arcpy.da.SearchCursor and arcpy.da.UpdateCursor to read and store a couple of values from Table 1 Field A, then insert those values as a single string within Table 2 Field B.  So basically 2 values into 1.  There's a lot more to it of course, but this is the main piece I'm struggling to figure out.  I think I need to use a dictionary to store the values, but I'm honestly not too familiar with them.  That along with some list comprehension.
import arcpy

# Local Variables:
Table_1 = "C:\Path\TestArea.gdb\Table_1"
Table_2 = "C:\Path\TestArea.gdb\Table_2"

with arcpy.da.SearchCursor(Table_1, "MILE_TEXT") as cursor:
    for row in cursor:
        #Do something to grab the values in the MILE_TEXT field and concatenate them in preparation for the UpdateCursor

with arcpy.da.UpdateCursor(Table_2, "MILE_MARKERS") as cursor:
    for row in cursor:
        #Do something here to take the values collected in the SearchCursor and update a single field/cell with that value
        cursor.updateRow(row)

print "Complete Without Errors"

Here's a highlight of the data I'm trying to capture and use to populate

Here's the single cell I'm trying to populate with the captured values


texts: Get correct measurements for fantasy world map

texts:
I read lots of questions about fantasy maps and custom crs but none of the answers could solve my problems.

Here is what I have:

A planet with an equatorial circumference of 38 211 km (radius: 6081.5 km)
A bunch of shape files for the continents and some islands as well as country borders and so on
A raster graphic in the background, already georeferenced
Shape files and graphics coordinates in longitude latitude degrees.

Here is what I want to get:

a world map in a cylindrical projection
a custom crs to get metric measurements from the degree coordinates of my layers.

Here is what I tried:
 - define custom crs:

+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6081500 +a=6.08146953 +b=6.08146953 +units=m +no_defs

(+a=r[km]/(2π*1000) and +b the same thing)

Then I clicked on project --> properties --> "CRS" and selected the new crs
Under "General" --> "Measurements" --> "Ellipsoid" I chose "Custom" and entered the radius of my planet in meters for both semi-major and semi-minor values.
I right-clicked the layers and chose the same crs for them as for the project.
I selected the "measure line"-tool and measured a line on the equator that should be around 38211 km or 360°. It never was. I measured a line on the poles that should be around 0 km, never got close to 0.

I changed the +a and +b values in the crs, I made sure for my shape files to have the right coordinates, but I can't find any errors.

So the final question is:
How would a crs look like, that achieves the following parameters for my custom planet map:

Cylindrical projection
Poles at +90° (Northpole) and -90° (Southpole), 9,552.798 km (quarter circumference of the sphere) away from the Equator
expanding -180° and +180° (119,105.596 km, half circumference) from west to east?


texts: Creating bathymetry surface from xyz in GRASS?

texts:
I am trying to create a 3D bathymetry surface from xyz data.

There is approximately 300,000 points with lat, long and z (depth data). The data is sonar point data similar to LiDAR data.

I want to use GRASS and I'm looking around for the right functions.

I can import the data to grass with v.in.ascii -z and I tried using v.extrude which killed my computer (must be the number of points).

I might be able to create contours in QGIS and use v.to.rast and r.surf.contour but I think I'll lose resolution if it works at all.

I'm not familiar with GRASS and I guess that, this is probably a relatively easy thing to do but I've hit a wall with it. I would imagine someone has encountered this before.


texts: Error in .BIOMOD_RangeSize.check.args(CurrentPred, FutureProj, SpChange.Save) : CurrentPred & FutureProj dimensions mismatched!

texts:
Colias_bin_proj_current <-
  stack(
    c(
      ca = "Colias.erate/proj_current/individual_projections/Colias.erate_EMcaByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img",
      wm = "Colias.erate/proj_current/individual_projections/Colias.erate_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
    )
  )

Colias_bin_proj_2050_CC45 <-
  stack(
    c(
      ca = "Colias.erate/proj_2050/individual_projections/Colias.erate_EMcaByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img",
      wm = "Colias.erate/proj_2050/individual_projections/Colias.erate_EMwmeanByTSS_mergedAlgo_mergedRun_mergedData_TSSbin.img"
    )
  )

crs(Colias_bin_proj_current) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
crs(Colias_bin_proj_2050_CC45) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

resample( Colias_bin_proj_current, Colias_bin_proj_2050_CC45, method='bilinear' )

SRC_current_2050_CC45 <-
  BIOMOD_RangeSize( Colias_bin_proj_current,
                    Colias_bin_proj_2050_CC45
                      )

  Colias_bin_proj_current
class      : RasterStack
dimensions : 329, 272, 89488, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.04166667, 0.04166667  (x, y)
extent     : 65.58333, 76.91667, 24.16667, 37.875  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
names      : ca, wm
min values :  0,  0
max values :  1,  1

 Colias_bin_proj_2050_CC45
class      : RasterStack
dimensions : 79, 101, 7979, 2  (nrow, ncol, ncell, nlayers)
resolution : 0.1666667, 0.1666667  (x, y)
extent     : 61, 77.83333, 23.83333, 37  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
names      : ca, wm
min values :  0,  0
max values :  0,  0

I have resample the current and 2050 but I got same dimension error and my dimensions are still different. How can I solve this issue?


texts:
Link to GEEScript

var L8 = L8.filterBounds(geometry).first()
Map.addLayer(L8,{max:3000},'Before Reprojection')
var proj = ee.Projection('EPSG:4326').getInfo()
var L8r = L8.reproject(proj.crs,proj.transform)
Map.addLayer(L8r,{max:3000},'After Reprojection')

A reprojection to EPSG:4326 would be benefictial. But after Reprojecting, the Image seems Oddly Broken into 3 Large Pixles.

Adding a Picture Manifactured from the Projection will result in much more then 3 Pixels an a usefull looking resolution :

var coordspic = ee.Image.pixelCoordinates(ee.Projection('EPSG:4326')).select('x')
Map.addLayer(coordspic,{max:2,min:-7},'Coordinates')

What am I missing here ?


texts: Grouping database results based on intersection using MySQL GIS?

texts:
I am using MySQL GIS (5.7.18).

I have a number of locations in the database. I also have a number of other polygons (not in the database).

I would like to perform a single query, intersecting all the places with all the polygons, grouping the places based on the polygon that they intersect with.

If this is possible, how can it be done?


texts: Transforming to GPS coordinates using projections with osgeo with Python

texts:
I am trying to transform a shapefile and its datum to GPS (EPSG Projection 4326 - WGS 84) but I get very different result on Google Earth, I write:

from osgeo import ogr, osr

file_name = 'POIS.shp'
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(file_name)
layer = dataset.GetLayer()
in_spacial_ref = layer.GetSpatialRef()
print (in_spacial_ref)
out_spatial_ref = osr.SpatialReference()
out_spatial_ref.ImportFromEPSG(4326)
coord_trans = osr.CoordinateTransformation(in_spacial_ref, out_spatial_ref)

x, y = 734150.2990000024, 6465754.138
x_gps, y_gps, _ = coord_trans.TransformPoint(x,y, 0)

# Google earth  3° 26.139'E, 45° 17.418'N,
print (x_gps, y_gps)
# 3.435642392921613 45.290299860385545

print (in_spacial_ref) gives me:

PROJCS["RGF_1993_Lambert_93",
    GEOGCS["GCS_RGF_1993",
        DATUM["Reseau_Geodesique_Francais_1993",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",700000.0],
    PARAMETER["False_Northing",6600000.0],
    PARAMETER["Central_Meridian",3.0],
    PARAMETER["Standard_Parallel_1",44.0],
    PARAMETER["Standard_Parallel_2",49.0],
    PARAMETER["Latitude_Of_Origin",46.5],
    UNIT["Meter",1.0]]

How can I get the exact values of Google Earth? What I am missing here?


texts:
I've created script in Python for ArcGIS to save shapefile as layer file... the code throws an error:

  ERROR 000732: Input Layer: Dataset D:\core i 5 data\symon assignment\Swat\clipped\buffer.shp does not exist or is not supported

Here is my code:

import arcpy
from arcpy import env

env.workspace = "D:/core i 5 data/symon assignment/Swat/clipped"
buf= "D:\core i 5 data\symon assignment\Swat\clipped\buffer.shp"
ou2= "D:\core i 5 data\symon assignment\Swat\clipped\layer1.lyr"
arcpy.SaveToLayerFile_management(buf, ou2, "ABSOLUTE")


texts: disable onEachFeature at zoom level Leaflet

texts:
In there a way to disable onEachFeature at zoom level? I have a GeoJSON object with highlightFeature and resetHighlight functions. I want to disable resetHighlight at zoom level 8

    var rejon_3;
    function resetHighlight3(e) {
        rejon_3.resetStyle(e.target);
        info.update();
    }

    function onEachFeature33(feature, layer) {
        layer.on({
            mouseover: highlightFeature
        });
    }

    function x3(e) {map.fitBounds(e.target.getBounds())}

    function onEachFeature3(feature, layer) {
        layer.on({
            mouseover: highlightFeature,
            mouseout: resetHighlight3,
            click: x3
        });
    }

    rejon_3 = L.geoJson(exp_3, {
        style: style3,
        onEachFeature: onEachFeature3
    }).addTo(map);

    feature_group.addLayer(rejon_3);


texts: Removing Interior and exterior polygons using arcpy

texts:
I have a polygon layer that has donuts and also exterior parts.
Using python, I am able to get access to each part, but  what I would like to know is using arcpy, once I have a part, how do I go about deleting all the vertices for that part.
This is what I have

and this is what I would like to end up with


texts: Cannot add capital N (letter N) to attribute field in QGIS?

texts:
New to QGIS and couldn't find anything related to this.

If I want to add a capital letter n, "N", to an attribute field QGIS does not accept it. The shapefile I am editing has 'Y' and 'N' in a text field. I have a feeling this is related to the NULL that is automatically applied to empty fields, but have no idea how to resolve.

The only way I've been able to add a capital N to a field is via the field calculator every time, which is too time consuming for the amount of editing I am doing.


texts:
With an enormous amount of mapfiles in one directory, is it possible to return a list of available mapfiles?

I run  MapServer in a docker container:

docker run -d -p 8181:80 -v $MAPDATA:/map bjoernschilberg/docker-alpine-mapserver

A single mapfile is available with:

localhost:8181/cgi-bin/mapserv/?map=/map/1990.map&service=WMS&version=1.3.0&request=GetCapabilities&


texts: Field added to raster data is not visible in attribute table

texts:
I am new to raster data and my background is software engineering so this might be a trivial question for GIS specialists.

Some raster files have attribute table but some do not.

On those with attribute tables, I can add a field in that table using Add Field (Data Management) tool but I cannot see the field in the attribute table.

In fact, when I identify any pixel by Identify tool, it does not show the field name so I am wondering where the field is gone?

EDIT :
I have tried to do Build Raster Attribute Table but it produces error. My raster data's pixel type is: signed integer. Steps I am doing now:

raster to polygon conversion(input:inpRaster, output:resultPolygon)
polygon to raster conversion(input: resultPolygon, output: resultRaster)
add field(input raster: resultRaster)
no error occurred. Shows field added message. But does not show it in the attribute table. If I remove the raster layer from arcMap and add it again from arcCatalog, then it shows the field.

When I try to create another new field with the same name, it says that field already exists. So why doesn't it show all fields in the attribute table?


texts: Creating Custom Symbology Patchs for ArcGIS Pro Map Legends?

texts:
Does anyone know if it is possible to create a custom symbology patch to be applied to a map legend in ArcGIS Pro?

I know this was possible in ArcGIS Desktop within the legend properties dialog, Items tab, style properties window.


texts:
I'm trying to read many ASTER images (1B mainly) in IDL (I can try in R also). I'm used to work with Landsat, but the aster metadata is complicated. I write a code with a lots of conditional that can read the image parameters (angles, dates, etc). But I want to know if there is other way I can read this.

Besides, I want to know if all the aster metadata have the same format, so I can batch all the images with the same code (or I need to modify it image by image).


texts: Convert Sentinel 2 1C product from reflectance to radiance

texts:
Does anyone know the equation for converting Sentinel 2 level 1c products from reflectance to TOA radiance?

I've looked at the Sentinel user guides but can't find anything. I'd rather use ArcGIS than SNAP, if possible.


texts: Openlayers 3 styling dynamic radius feature performance issue

texts:
I am migrating from 2.14 and have some features that come with a radius property.

The expected behaviour is that when the users zooms in or out the feature renders without perceptible jumps as with 2.14.

http://gyazo.com/bb46051994981d6d22cec4f7f88ff03a (ol 3)
http://gyazo.com/4f6a2793abacc597029c3c28cfec316f (ol 2.14)

Current code (v3)

var layer = new ol.layer.Vector({
    source: new ol.source.Vector(),
    style: styleFunction
});

var styleFunction = function(feature, resolution) {
    return [new ol.style.Style({
        stroke: new ol.style.Stroke({
          color: 'red',
          width: 3
        }),
        fill: new ol.style.Fill({
          color: 'rgba(255, 0, 0, 0.1)'
        }),
        image: new ol.style.Circle({
            radius: calculateRadius(feature, resolution),
            fill: new ol.style.Fill({
              color: 'rgba(255, 0, 0, 0.1)'
            }),
            stroke: new ol.style.Stroke({
                color: 'orange',
                width: 2
            })
          })
      })];
};

Previous code (v2.14)

var context = {
    getRadius: function(feature) {
        var radius = setup.DEFAULT_RADIUS;
        if (feature.attributes && feature.attributes.radius) {
            radius = feature.attributes.radius;
        }
        return radius / feature.layer.map.getResolution();
    }
};

var defaultStyle = new OpenLayers.Style({
  "pointRadius": "${getRadius}"
}, { context: context });

var layer = new OpenLayers.Layer.Vector("layer", {
    styleMap: defaultStyle
});


texts: SQL / Select by Attributes query for "does not include text found in another specified field"

texts:
I'm using ArcMap 10.2 and I'd like to write a query to select features where Field A does not include text found in Field B.

Specifically, Field A contains street numbers and street names. Field B contains only street names. I'd like to select features where Field A does not include the street names in Field B.

So far I've tried -

"Field_A" NOT LIKE'% "Field_B"%'


texts:
I would like to combine several catchment areas (points feature) into one catchment area. There are examples of my catchment areas.

Each area has one source (POI).I will use the combined area to create a map like this one(isochrone map by Anita Graser)

What should I  do to get my areas combined?


texts: Uploading 3D Buildings with Textures to ArcGIS Online

texts:
I'm doing a small project where I want to create 3D buildings from building footprints. I have acquired the heights of each building and have extruded them in ArcScene. Then I use Layer 3D to feature class to create a multi patch then exported it to COLLADA using the tool multi patch to COLLADA.

Then I took into SketchUp and have colorized the buildings using its photos taken using a camera. Now I want to use some software or something open source to view my buildings. I was told ArcGIS online can do that. I was searching and I couldn't find any option where I can import the texture COLLADA. There were only options like CSV, zip files.

Any suggestions?


texts:
I'm fairly new to using Survey123 and have ran into a problem. I've created a multiple choice question where you can select more than one answer. With each answer, I would like to show a specific relevant sub-question. Basically, if the answer is A, I want it to show only A sub question. If it's A and C, only show A and C sub questions. If it's C, only show C... So far, I have used the "${question1} = '1' or ${question1} = '2'" expression but that doesn't work.

How can I accomplish this?


texts: ESRI ArcGIS routing API extremely slow compared to other providers

texts:
I'm implementing a wrapper that requests route information from different geospatial / GIS API providers. Input is always two or more waypoints (latitude, longitude pairs).

It is not my goal to do a public benchmark or a comparison between different service providers. The problem is, it works pretty well with Google Maps, Bing maps and Here maps but requesting a route from Esri route API takes forever.

For all providers, I tried to optimize settings as far as I could:
- Limiting the requested output to geometry, time and distance
- Using just default driving modes (such as car, no real time information...)

EDIT:

Request URLs and response times (requested from Berlin, Europe) below. I executed the requests multiple times, even at different days and the times are almost the same. Coordinates are always Paris to Prague.

Esri: 46,26 seconds(!)

https://route.arcgis.com/arcgis/rest/services/World/Route/NAServer/Route_World/solve?f=json&returnDirections=false&token=<my_api_token>&stops=2.3522,48.8566;14.4378,50.0755

Bing: 611 milliseconds

https://dev.virtualearth.net/REST/V1/Routes/Driving?key=<my_api_token>&ra=excludeItinerary,routePath&wp.1=48.8566,2.3522&wp.2=50.0755,14.4378

HERE: 664 milliseconds

http://route.api.here.com/routing/7.2/calculateroute.json?xnlp=CL_JSMv3.0.17.0&app_id=<my_api_token>&mode=fastest%3Bcar&representation=display&waypoint0=geo!48.8566%2C2.3522&waypoint1=geo!50.0755%2C14.4378&jsoncallback=H.service.jsonp.handleResponse(5)

Google: 92 milliseconds

https://maps.googleapis.com/maps/api/js/DirectionsService.Route?5m4&1m3&1m2&1d48.8566&2d2.352200000000039&5m4&1m3&1m2&1d50.0755&2d14.437800000000038&6e0&12sen-US&23e1&callback=_xdc_._95pqrj&key=<my_api_token>
&token=117312

Has anyone experienced the same? Do I miss any steps?


texts:
One of the columns in my attribute table is called "table_ot".
Some of the features in that column start with a letter. Example: table_ot = G 190
I want to delete all the letters at the start of the feature's value and end up with something like this:
table_ot = 190

Can you please write down the correct code for the field calculator? I tried some things that don't run.


texts: Use QGIS plugins in custom standalone application

texts:
I am trying to create a QGIS custom standalone application. I want to add the Google Physical map from OpenLayers plugin to my application when it starts, but I cannot seem to be able to access qgis.utils. These are my imports:

import qgis.core
import qgis.gui
import qgis.utils
import PyQt4.QtCore
import PyQt4.QtGui

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.utils import *
import sys
import os

But when I try to print qgis.utils.plugins I get an empty dictionary; {}

In the Python console in QGIS, this code does what I want, I just can't seem to get it working in my custom app;

olplugin = qgis.utils.plugins['openlayers_plugin']
ol_gphyslayertype = olplugin._olLayerTypeRegistry.getById(0)
olplugin.addLayer(ol_gphyslayertype)


texts: USGS Relief map (GeoTIFF raster) misaligns with vector layers (QGIS)

texts:
I'm trying to correct the alignment of three layers (described below), but I'm not entirely sure how to go about re-projecting raster layers to my other layers. Ultimately my goal is to clip the relief map to the park boundaries to look similar to Grossenbacher's map. There are some park boundaries to remove/exclude, but I first want to make sure everything is aligned.

The main thing that is confusing me is the many different projections of my layers, and which one to choose as the main projection. As well as correctly re-projecting the raster layer.

Layer 1 (brown, TIGER census tracts)

CRS: EPSG:4269 - NAD83 - Geographic
Source: https://www2.census.gov/geo/tiger/TGRGDB16/

Layer 2 (yellow, National Park boundaries)

CRS: EPSG:3857 - WGS 84 / Pseudo-Mercator - Projected
Source: https://public-nps.opendata.arcgis.com/datasets/national-park-service-park-unit-boundaries

Layer 3 (gray, USGS Relief map [Shaded Relief Land - Gray - Conterminous United States 200 Meter Resolution])

CRS: Generated CRS (+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs) - Projected
Source: https://nationalmap.gov/small_scale/atlasftp.html#srgy48i


texts: ogr2ogr error loading KML file in PostGIS with PHP

texts:
So I have a KML file that I want to save to my PostGIS database,
Currently I have the following code using PHP:
$file = $request->file("item");

$name   = time();
$name   = $name . '.' . $file->getClientOriginalExtension();
$temp = 'uploads/' . $user_app->id . '/kml/';
$path = public_path('uploads/' . $user_app->id . '/kml/');
if (!file_exists($path)) {
    mkdir($path, 0777, true);
}
$path   = public_path() . '/uploads/' . $user_app->id . '/kml/';
$command = "ogr2ogr  -f PostgreSQL PG:\"host=".$HOST." dbname=".$DBNAME." port=".$PORT." user=".$USER." password=".$PASSWORD."\"   ".$path.$name."  -nln siap 2>&1";

Log::debug($command);
$res = shell_exec($command) ;
Log::debug($res);

This is the command line,
ogr2ogr  -f PostgreSQL PG:"host=127.0.0.1 dbname=postgis_24_sample port=5432 user=postgres password=123456"   C:\xampp\htdocs\siap-padron\public/uploads/2044492/kml/1597284059.kml  -nln siap 2>&1

If I run it in PowerShell or in the system symbol it inserts it correctly in the PostGIS base, but when it is run in PHP it gives me the following error:
local.DEBUG: FAILURE: Unable to open datasource `C:/xampp/htdocs/siap-padron/public/uploads/2044492/kml/1597284059.kml' with the following drivers.
-> `JP2ECW'......

Any comment would be very helpful.


texts: How to make ui.Panel elements transparent (Earth Engine)

texts:
I would like to be able to adjust the transparency of the Earth Engine ui.Panel element. The documentation for panel style has no apparent option to set opacity/transparency; is there any way to do this?

Here is a panel example:

var panel = ui.Panel({
  style: {
    height: '235px',
    width: '600px',
    position: 'bottom-right',
  }
});

Map.add(panel);


texts: Quick OSM error "not attribute 'get_bounding_box'"

texts:
I'm new on Qgis.
I want to make query using Quick OSM pluging in QGIS 3.16. When I generate a query, it's returning me this error : AttributeError: 'Dialog' object has no attribute 'get_bounding_box'
But when I execute the query, it's working
Any idea what is happening ?
query :
[timeout:2500]
[maxsize:2000000000];
area[admin_level='6'][ref='90']->.area;
rel(pivot.area)->.rel;
(
node[highway](area.area);
node[highway](r.rel);
way[highway](area.area);
way[highway](r.rel);
rel[highway](area.area);
rel[highway](r.rel);
);
(._;>;);
out;


texts: Problem with grouping items in print composer- QGIS 2.14

texts:
This is print composer:

when i group all elements with "group items" and minimize the frame- some items step out from the frame:

How can i group all items, and still maintain the same ratio of distances between all items, so all items remain in their position?

In ArcGIS it is possible.


texts: High res monthly satellite maps including ice cover

texts:
I have some old month by month NASA satellite maps showing the world in 2004 at 5400x2700px. I've noticed NASA has newer cleaner versions of these maps at http://neo.sci.gsfc.nasa.gov/view.php?datasetId=BlueMarbleNG-TB&date=2004-12-01, however the maximum resolution is only 3600x1800. Also both map sets have ice cover removed.

Are there any month by month raster satellite world map sets out there with a resolution of 5400x2700px or higher and including ice cover?


texts: Troubleshooting batch reprojection script that fails some of the time without error but works other times?

texts:
This script, up until now, has successfully batch reprojected datasets of rasters to designated projections. The results of running my script on my new dataset is that the script successfully reproject the first 6 out of 200 and then finishes without an error. How is the best way to troubleshoot this issue?

import arcpy, sys, os

InFolder = r"C:\TX_DataSets\Harris 1944\Copy"
OutFolder = r"C:\TX_DataSets\Harris 1944\Reprojected"
OutSR = arcpy.SpatialReference(26915) # UTM 15

arcpy.env.workspace = InFolder
print "Starting reprojecting script"

rasterList = arcpy.ListRasters("*","tif")
for Ras in rasterList :
     inRaster = os.path.join(InFolder, Ras)
     spatial_ref = arcpy.Describe(inRaster).SpatialReference
     if spatial_ref.type <> "Unknown":
          outFileName = "{}.{}".format(os.path.splitext(Ras)[0], "tif") # Remove .sid and add .tif to Ras filename
          outRaster = os.path.join(OutFolder, outFileName) # Path and Filename of output raster (.tif)
          arcpy.ProjectRaster_management (inRaster, outRaster, OutSR)
          print "reprojecting " + str(outRaster)
print "Finished Reprojecting"


texts: How to cluster interconnected lines in Python without PostGIS?

texts:
I'm working with a subset of data from OpenStreetMap, where I have these island-like groups of line features. The groups are not connected to one another, but the lines are all interconnected within each group. I'd like to aggregate groups of line features that 1) are contiguous and 2) share a common attribute value, into single features, with non-grouping attributes aggregated.

This may sound like a dissolve operation, but a dissolve will typically aggregate all features in the dataset with a common value into a single multipart feature, irrespective of spatial relationship. Rather, I want only features that are connected to each other (directly or by intermediate features) to be aggregated into an output feature.

PostGIS's cluster functions offer exactly what I want, except for the database dependency. In PostGIS, I would do a query like this:

SELECT ST_Collect(the_geom) AS geom,
       min(osmid::integer) AS osmid,
       string_agg(distinct highway, ',') AS highway,
       string_agg(distinct ref, ',') AS ref
FROM
    (SELECT *,
            ST_ClusterDBSCAN(the_geom, 0, 1) OVER() AS _clst
     FROM mytable) q
GROUP BY _clst
ORDER BY osmid

Is it possible get this same result without PostGIS, using only Python libraries like Shapely, GeoPandas, or similar?


texts: How to get X and Y coordinate values of the map?

texts:
I want the upper right and lower left coordinate of the exposed portion of the map, and for the whole map.

What is the way to do this?


texts: using ogr2ogr to import a shapefile into a PostgreSQL/PostGIS database

texts:
I am trying to load a shapefile into a PostgreSQL/PostGIS databases with the following statement run within the OSGeo4W shell
ogr2ogr -f "PostgreSQL" PG:"host=localhost port=5432 dbname=testing user=postgres password=12345" "D:/TEST_RUBBISH/mains_topo.shp"

I get the following warning which effectively preempts the failure of the insert operation

Warning 1: Geometry to be inserted is of type Multi Line String, whereas the layer geometry type is Line String.
Insertion is likely to fail

I have been all over the internet and into the 'PostGIS in Action' book to try and work out how to resolve this but no joy so far
Could anyone out there suggest a solution


texts:
I am using ArcGIS Desktop 10.4.

I am trying to iterate a process through a Feature Class to calculate individual (and unique) values in a particular field. I can calculate the field just fine, but it ends up recalculating the new unique value for all of the features in the feature class instead of each individual row. I cannot convert to feature layer or layer file because I need to add the data to the current Feature Classes without adding more Feature Class to that geodatabase. (I am already iterating it through 231 Feature Classes and do not need more...) I have tried Search Cursor already, and could not find any other code because arcpy.IterateFeatureSelection() is meant for model builder and not python.
I was expecting my update cursor in "fc_update in fc_updates" would provide that service, but it did not. Otherwise it prints and defines the data exactly how I intent to use it.

I left out the rest of my code as it is not pertinent to the part that needs to select the data.

My code is below:

fc_updates = gp.UpdateCursor(fc)

for fc in fcs:
    if fc.startswith("GP_"):
        print "Start Feature Class at: " + str(datetime.datetime.now().time().strftime('%H:%M:%S'))
        print "Feature Class: " + fc
        #arcpy.AddField_management(fc, field_name,"Short")

        for fc_update in fc_updates:
            obj_id = str(fc_update.OBJECTID)
            flt_desig = str(fc_update.flt_line)
            imag_start = str(fc_update.Start_Imag)
            imag_end = str(fc_update.End_Imag)
            n = n + 1
            print "-"*30
            print "Attribute Information . . ."
            print "image start = " + str(imag_start)
            print "image end = " + str(imag_end)
            imag_entry = (str(imag_start) + "," + str(imag_end)).split(",")
            print "image entry = " + str(imag_entry)
            start = int(imag_entry[0])
            end = int(imag_entry[1])
            express1 = start-n
            express2 = start+n
            print (str(flt_desig), str(imag_start), str(imag_end))
            print "-"*30


texts:
QGIS-Version 2.14.5 LTR, QGIS Server running on Windows 2008 R2

I encounter a '500 Internal Server error' when using the url http://localhost/qgis-ltr/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities&map=flst_ster.qgs

My httpd_qgis-ltr.conf looks like this:

LoadModule fcgid_module modules/mod_fcgid.so

SetEnv PATH "C:\OSGeo4W\\bin;C:\OSGeo4W\\apps\qgis-ltr\bin;C:\OSGeo4W\\apps\grass\grass-@grassversion@\bin;C:\OSGeo4W\\apps\grass\grass-@grassversion@\lib;C:\Windows\system32;C:\Windows;C:\Windows\System32\Wbem"
SetEnv QGIS_PREFIX_PATH "C:\OSGeo4W\\apps\qgis-ltr"
SetEnv QT_PLUGIN_PATH "C:\OSGeo4W\\apps\qgis-ltr\qtplugins;C:\OSGeo4W\\apps\qt4\plugins"
SetEnv TEMP "C:\Users\sjo\AppData\Local\Temp"
SetEnv PYTHONHOME "C:\OSGeo4W\apps\Python27"
SetEnv PYTHONPATH "C:\OSGeo4W\apps\qgis-ltr\python"

Alias /qgis-ltr/ C:\OSGeo4W\/apps/qgis-ltr/bin/

<Directory "C:\OSGeo4W\/apps/qgis-ltr/bin/">
    SetHandler fcgid-script
    Options ExecCGI
    Order allow,deny
    Allow from all
</Directory>

I tried to replace (as suggested in various sources) to replace fcgid by cgi but this causes the same error and an additional Premature end of script headers: qgis_mapserv.fcgi.exe in the apache error.log

I remember getting this to work with 2.8.9.LTR (here I HAD TO replace fcgid by cgi) so have I missed any changes? What might have gone wrong here?


texts: How to make Hierarchical Tree diagram from polygon layer?

texts:
I have a polygon layer with 5 hierarchical levels. I added 3 fields to indicate the level and the group. The hierarchy is also shown spatially, ie. a big polygon contains the other lower level polygons that are contained in that group.

What is the best tool to create diagram of hierarchical tree from the data?

I don't have license for Schematics or Visio.  I have Excel 2007.


texts: Creating closest line along the point''s azimuth using PostgreSQL

texts:
I have a point table which contains start point and end point:

When I am writing this query in PostgreSQL to create line it gives me this output.

select distinct on (a.id) a.id,
        min(st_astext(st_makeline(a.geometry, b.geometry))),
        min(st_distance(a.geometry, b.geometry)) as dist
into results
from fpts1 as a,
     fpts1 as b
where not a.geometry = b.geometry
      and not a.link_id = b.link_id
      and ST_DWithin(a.geometry, b.geometry,0.0002)
group BY a.id,
         ST_Distance(a.geometry, b.geometry);

How I can avoid line creation like this and get output like below image.


texts: Google Maps Android SDK - Get features based on user location

texts:
I am using Google Maps Android Utility SDK, using GeoJSON features.

I would like to, based on the user Location (which I already have) determine the features that contain this point.

Does Google Maps provide such utility or should I implement by myself?


texts:
how can I add OSM be the base layer? Below is what I have and the part within asterisk is where I want to place an OSM.

var mapPanel = new GeoExt.MapPanel({
region: "center",
center: [27.795971, -25.609990],
zoom: 11,
layers: [ * * new OpenLayers.Layer.WMS("OpenLayers WMS",
        "http://labs.metacarta.com/wms/vmap0", {
            layers: 'basic'
        }, * * {
            buffer: 0,
            // exclude this layer from layer container nodes
            displayInLayerSwitcher: false
        }),
    new OpenLayers.Layer.WMS("Madibeng City",
        "http://localhost:8080/geoserver/cite/wms", {
            layers: "topp:Madibeng0",
            transparent: true,
            format: "image/gif"
        }, {
            buffer: 0
        }),
    new OpenLayers.Layer.WMS("Madibeng Road Network",
        "http://localhost:8080/geoserver/cite/wms", {
            layers: "Madibeng_LM_RoadNetwork0",
            transparent: true,
            format: "image/gif"
        }, {
            buffer: 0
        }),
    new OpenLayers.Layer.WMS("Madibeng Paved Locations",
        "http://localhost:8080/geoserver/cite/wms", {
            layers: "topp:Verf_Paved0",
            transparent: true,
            format: "image/gif"
        }, {
            buffer: 0
        })
 ]
 });


texts: QGIS Server does not publish Oracle Spatial Layer

texts:
I want to publish a Oracle Spatial Layer which is defined in a QGIS Project with the QGIS Server as a WMS.
The capabilities document of the QGIS Server does not include the layer description of the Oracle Spatial Layer. If i try this for example with a PostGis Layer it will work. Is there a solution about this.

My Version of QGIS is 2.4 for Windows 64 bit and i installed the QGIS Server over the OSGEO setup on the same system.


texts: How do I get polygon shapefiles to project correctly in WGS 84?

texts:
I have a map with three layers:

"Dorset Points", point shapefile
"Counties...", polygon shapefile
"Kilwood Park", raster

I'm trying to project all three layers onto WGS 84. The project CRS is set as WGS 84 and I've manually set each layer to project to WGS 84. I've enabled "on the fly projection" as well. However, layers 2 and 3 align properly while layer 1 lies out in the ocean. Also, when scrolling with my mouse, my coordinates are displayed as what looks to be northings and eastings instead of decimal degrees as I would expect with WGS 84. When I zoom to layer 1, my scrolling mouse shows decimal degrees. I know layer 2 has decimal degrees as you can see in the attribute table. I've tried using the field calculator to produce x and y coordinates using $x and $y but my output was null.

Please let me know of anything else to try! I know there are many resources on projecting but so far I have tried every fix I could find.


texts: Configure new SQL view doesn't appear for different DB's in the same PostGIS Database

texts: I have PostGIS Database with different DB's. Each one of them has different workspace:store in Geoserver. When I'm opening Create new layer page and picking different workspace:store for some of them option "Configure new SQL view " appears and for some of them not. I could not find any correlations DB/workspace/store definition and this behavior... Any thoughts?