texts: PostgreSQL+PostGIS creating GeoJSON with ST_MakeLine

texts:
Trying to get a valid geojson response from my PostgreSQL(9.6.3) + PostGIS (2.3.2) to track routes. I have a table of GPS positions of the form:

CREATE TABLE public.gps
(
  id integer NOT NULL DEFAULT nextval('gps_id_seq'::regclass),
  device_id integer, -- Device ID of GPS Device
  timecol bigint, -- Timestamp we use for linking video and GPS
  distance_from_hub double precision, -- Distance from nearest hub
  distance_from_hub_id integer, -- ID of nearest hub
  itudes_3d geography(PointM,4326)
)

I am trying to get a valid GeoJSON with the route as a LineString. The coordinates of each point are stored in the itudes_3d. I am querying with:

select jsonb_build_object(
        'type', 'FeatureCollection',
        'features', jsonb_agg(feature)) AS route
        FROM ( select jsonb_build_object(
            'type','Feature',
            'properties',jsonb_build_object(
              'opa', 1,
              'fillopa', 0.5,
              'type', 0,
              'weight', 1,
              'start_time', :start_time:,
              'end_time', :end_time:,
              'colour',
              '#333'),
            'geometry', geometry.geoline::json)
            FROM (select ST_AsGeoJSON(ST_MakeLine(geom.itudes_3d))::jsonb as geoline
              FROM (select itudes_3d::geometry, ST_M(itudes_3d::geometry) as gtime
                FROM gps
                WHERE device_id=:device_id: AND
                timecol > :start_time: AND
                timecol < :end_time:
                ORDER BY timecol ASC)
              AS geom)
            AS geometry)
         AS feature

where :start_time:,:end_time: and ":device_id:" are placeholders for this explanation.

What I get back is:

{"type": "FeatureCollection", "features": [{"jsonb_build_object": {"type": "Feature", "geometry": {"type": "LineString", "coordinates": [[174.821949, -35.829088], [174.820518, -35.82855], [174.818997, -35.827884], [174.818997, -35.827884], [174.818696, -35.827857], [174.818363, -35.827694], [174.816604, -35.826887], [174.816275, -35.826737], [174.815868, -35.826569], [174.815188, -35.826255], etc...

The problem I have is the {"jsonb_build_object": label in the features... really I want:

{"type": "FeatureCollection", "features": [{"type": "Feature", "geometry": {"type": "LineString", "coordinates": [[174.821949, -35.829088], [174.820518, -35.82855], [174.818997, -35.827884], [174.818997, -35.827884], [174.818696, -35.827857], [174.818363, -35.827694], [174.816604, -35.826887], [174.816275, -35.826737], [174.815868, -35.826569], [174.815188, -35.826255], etc...

Can anyone identify where I have gone wrong in constructing this?

To me it looks like the subquery is being named in the output.


texts: How to count point shp inside polygon shp using Python?

texts:
This is my first time here, so I'll try to be as specific as I can.

What I'm trying to do is found out which polygon contains the most points and write the result in a text file.

Here is my script : (tab25Select are the polygons and chefLyr are the points)

chefLyr = arcpy.MakeFeatureLayer_management(cheflieuIN, "chefLyr")

with arcpy.da.SearchCursor(tab25Select, ["Numero", "Name"]) as rows:
    for row in rows:
        chefTab25Max = 0
        chefTab25 = arcpy.GetCount_management(arcpy.SelectLayerByLocation_management(chefLyr, "WITHIN", tab25Select))
        arcpy.AddMessage(cheflTab25)
        for chf in cheflTab25:
            if chf > chefTab25Max:
                chefTab25Max = cheflTab25
                a = row[0]
                b = row[1]

ficTxt.write("\r\n The polygon containing the most point is :{0}".format(chefTab25Max))

ficTxt.close()

If anyone can help :)


texts: Using cursor to populate a new field by concatenating existing fields

texts:
I have a question for a current assignment that requires me to write a script that will add a new field to a feature class and then use a cursor to populate the new field by concatenating five other fields (name, address, city, stateabb, zipCode). The new field is FullAdd and will be the full mailing address in one field from the information that is already in the table.

I need help with the update cursor part of my code.

How do I make it write the concatenation to the new field?

I've tried so many things after the 'with cursor' statement and for loop and have no idea where to go next!

import arcpy
from arcpy import env
env.workspace = r"C:\mygdb.gdb"
env.overwriteOutput = True
# Set up search cursor for hospitals FC
# Make the hospitals FC a lyr
arcpy.MakeFeatureLayer_management("Hospitals", "Hospitals_lyr")

# Add the new field
arcpy.AddField_management("Hospitals_lyr", "FullAdd", "TEXT")

fields = ['NAME', 'ADDRESS', 'CITY', 'ST_ABREV', 'ZIP_CODE', "FullAdd"]

cursor = arcpy.da.UpdateCursor("Hospitals_lyr", fields)

try:
  with cursor as updateCursor:
    for row in updateCursor:
      name = row[0]
      address = row[1]
      city = row[2]
      stateAbb = row[3]
      zipCode = row[4]
      FullAdd = row[5]
      FullAdd = str(name + " " + address + " " + city + ", " \
      stateAbb     + " " + zipCode)
      updateCursor.updateRow(row)

  # Copy feature lyr back to disk
  arcpy.CopyFeatures_management("Hospitals_lyr", "Hospitals_2")
  print "done"

except:
    print "Error"


texts: I stored a googlemaps.tpk file in my sdcard. I am trying to run it with ArcGIS but it gives the error file is not a tiled package. My question is that how to run googlemap.tpk file with ArcGIS.


texts: Customizing QGIS print composer to use default layout created and custome table

texts:
I am doing a cadastral application as a plugin and I want to export some table in data [in a model] with the parcels displayed [ the geometry ] to pdf file.

And I am thinking to do that using  a QGIS print composer.But while creating layout

I was able to add loaded canvas
I was able to add attribute table but here I want my custom data table to be rendered
And I have added some lables and legened

Now what I really need is getting the loaded canvas with my data table in the layout prepared while clicking print composer.
what is the best way to get custome data to the print composer
I am able to code in QGIS iface enviroment and am looking solution in that aspect.


texts:
I've been having issues visualising data from a shapefile in GeoDjango using OpenStreetMap. I believe it's some sort of a scale issue because of how it appears (as you can see in the screenshot of the admin interface).

I've stored the data in a PostGIS database and set up the model using django as advised in the documentation. Another interesting thing to note is that loading up the data using QGIS works perfectly and looks exactly like the shapefile. Any ideas as to what I should try next?

I'd also be happy to post any code if people need to see it but I am new to GIS stuff and don't quite know what anyone would need to see right away.


texts:
I have two shapefiles:
1) one with several points with a buffer around each separated point
2) another which is a single polygon with several characteristics (therefore it has several polygons inside, with different sizes).

What I want to do:
I want to perform a "clip" and obtain the characteristics of the single polygon (2) within each point in (1).

This is what it looks like before the "clip":

And this is what is looks like after the "clip":

What I get is all the points together, but what I wanted was for them to be individual as well. How do I do that? Does anyone has any idea?


texts: Arcgis not seeing my column header to join tables

texts: I am trying to join a data table to the attribute table of a shapefile. I have done all the necessary steps to makes  sure it works: removed spaces, not use numbers in column headers, remove illegal characters etc. BUT still when I try to join the table, Arcmap is not reading the GEOID field that I am trying to use to join the two tables. I cannot see it listed in the fields to join. What could be the problem?


texts:
Does anyoe have a link to or know of a list of raster formats the above tool can work with? I'm on QGIS 2.8 using the latest version of the plugin.

I've tried a couple of different raster types and but have so far only had success on an .adf (ESRI grid) raster.


texts:
Can you provide me a little know-how of how you could code your streets as shown in the figure.
Initial class is assigned as A,B,C... and when certain street(say A) branches, it is coded , as A1, A2 with even and odd numbers corresponding to side the branching occurs. Similarly the number is added -either 1 or 2 as branching occurs. I believe network analysis with hierarcy is a good starting point.

I am okay with doing some Python and ArcPy.


texts: PostGIS ST_Covers doesn't match polygon after ST_Union

texts:
I have a simple table:
create table test (
    id serial,
    polygon geometry(Polygon, 3857)
);

And I want to insert two polygons in this table which should be created as points expand with intersection area between them:
insert into test(polygon)
values(ST_Transform(ST_Expand(ST_Transform(ST_SetSRID(ST_MakePoint(30.2833756, 50.4419441), 4326), 31997), 100), 3857));
insert into test(polygon)
values(ST_Transform(ST_Expand(ST_Transform(ST_SetSRID(ST_MakePoint(30.2841370, 50.4419441), 4326), 31997), 100), 3857));

After that I have a query to check if point is covered by polygons union:
with point_polygon as (
    select ST_Transform(ST_Expand(ST_Transform(ST_SetSRID(ST_MakePoint(30.2841370, 50.4419441), 4326), 31997), 100), 3857)
)
select ST_Covers(ST_Union(polygon), (select * from point_polygon))
from test
where ST_Intersects(polygon, (select * from point_polygon));

I can't understand why above query returns true for the first point and false for the second...
If I decrease the expand radius in select from 100 to 99.999 it returns true:
with point_polygon as (
        select ST_Transform(ST_Expand(ST_Transform(ST_SetSRID(ST_MakePoint(30.2841370, 50.4419441), 4326), 31997), 99.999), 3857)
    )
    select ST_Covers(ST_Union(polygon), (select * from point_polygon))
    from test
    where ST_Intersects(polygon, (select * from point_polygon));

Is it a rounding issue?


texts:
Alabama's data portal (portal.gsa.state.al.us) is down and I haven't come across any other GIS clearinghouses for AL that are as good as say, Florida's FGDL.

Does anyone know of a decent GIS clearinghouse specific to Alabama?


texts: Get neighbor census tracts for a given census tract

texts:
I have IDs for census tracts (USA census tracts). I want to get neighbor tracts for a given census tract through an API call.

There are APIs for get the census tract for a given location point (longitude, latitude). Like that, what are the available free APIs for getting the neighbors?

Example:

http://www.usboundary.com/Tools/Neighbors/Census%20Tract/California/San%20Francisco%20County/Census%20Tract%20260.02/440627


texts: How do I get the coordinates of the nearest point on a polygon?

texts:
The following SQL uses ST_Geometry to return a resultset of 10 records, ordered by distance from a fixed location (a point geometry) to the nearest point of the polygon. How do I get the coordinates of the nearest point on the polygon? For example, if a record shows up in the resultset and the distance is 100, I want to see the coordinates of the point 100 meters from my fixed location.

SELECT *
FROM
  (
    SELECT areas.XYZ, sde.st_distance(sites.shape,areas.shape)
    FROM ABC.site sites, XYZ.identifiedarea areas
    WHERE sites.id = 123
    AND sde.st_intersects(areas.shape, sde.st_buffer(sites.shape, 1000)) = 1
    ORDER BY sde.st_distance(sites.shape,areas.shape)
  )
WHERE ROWNUM <= 10 ;


texts: Is it possible from inside QGIS to analyze a number of occurrances (non-points) of a given variable repeating in a database to map it?

texts:
I have a database with each and every death in a country, each row in the database has a GEO_REFERENCE tag constructed by the concatenation of a code given to every division and subdivision instead of a GPS point, looking as follows:

 STATE_NAME    |  STATE_CODE  |    MUN_NAME    |  MUN_NAME  |  GEO_REF
Name of State  |      01      |    City Name   |   0001     |  010001

Already having a layer with the divisions and subdivisions named according to that GEO_REF code, Can it be counted how many times there has been a death at every one of this GEO_REF, in order to map it? Or should the database analysis must be done prior in a different program?


texts: Does a map service need write permissions for a file GDB?

texts:
If an ArcGIS Map Service is published to use data located in a file geodatabase on a network drive, it obviously must be able to read from those files to present the data to end users. However, does it also need write access to the file GDB? I could imagine this being required for creating or modifying lock files, for instance.

I've found two pages of ESRI's documentation for ArcMap, but it's not clear whether they apply to ArcGIS Server:

File geodatabases and locking doesn't explicitly say. The note about creating lock files "on access" would suggest perhaps it is needed, but this still isn't explicit.
File geodatabases and Windows Explorer says that you can set the file permissions to read only, but this is for ArcMap.

I'm using ArcGIS 10.3.

Side note

I'm well aware that serving up map service data from the network drive is typically not the best idea. Unfortunately, I am not the developer who designed or built the system, and I am tasked with deploying it without being given the opportunity to change this.

For future readers: If you have a choice, don't store map service data in a file GDB on a network drive. It's not good for performance.


texts: Get the N'th feature from QgsFeatureIterator fast

texts:
I have a vector layer with say 100.000 features. I would like to get the N'th feature from the layer when sorted by an attribute field name:

Feature 1
Feature 2
    :
Feature N   <--
    :
Feature 99.999
Feature 100.000

Using a QgsFeatureRequest with layer.getFeatures(...) I can sort the features:

request = QgsFeatureRequest()
request.addOrderBy(field_name)

feature_iterator = layer.getFeatures(request)

Then I can extract the N'th feature by:

feature = list(feature_iterator)[N - 1]

It takes some time to make a list of 100.000 features. I can limit the number of features since I don't need to request the features coming after N:

request.setLimit(N)

But if I need feature N=99.999 I still need to request almost all of the features.

Question:
Can I get the N'th feature more efficiently?


texts:
How could I perform a Tabulate Area using GDAL and Numpy?  I am trying to summarize how many pixels of each class in a Landcover raster are within each polygon of a Census Blocks shapefile.

I've tried playing with numpy.histogram2d, but didn't have any luck.


texts: Rasters not displaying correctly when opening from a project in QGIS

texts:
I have a few rasters of Landsat 8 imagery I made in GRASS and exported as geotiff files. I can load them fine in QGIS and they display correctly.

But if I save it as a project, then close and reopen QGIS and load the project, the rasters' colour maps are all wrong. There will be lots of white pixels and pixels with wrong colours. Sometimes QGIS crashes randomly with the project open. Removing the layer and re-opening them fixes it. I also notice when I go properties for one of these layers, it takes about 10-15 seconds for the window to come up.

Using 2.16.3. Rasters made in GRASS 7.0.4, with no compression, UInt16 and GDAL colour table. I think it's the colour table that is causing the problem, but I want the colour scheme preserved in the file.


texts: Alternative satellite and hybrid layers for Europe published on WMS

texts:
For an application I am building I require to show satellite (high resolution) and street view layers, preferably also a hybrid view. Unfortunately, due to the application not running on a typical GUI, restrictions of the Google Terms prohibit me from using Google Maps.

Looking for alternatives, I have had very little success in finding public WMS hosting alternative datasets. Hosting the data ourselves would require a redefinition of the project's scope which is, at this point, not an option.

Looking at questions asking for lists of WMS, I have unfortunately found many servers that are no longer online or providing the services I need. Do you know any WMS hosting the layers I need or, in case you don't, what WMS list would you consider most up to data and regularly maintained?


texts: How do I manually classify graduated data in QGIS?

texts:
Is there any way to manually classify Graduated Data in Qgis. In ArcGIS I always used to apply first the natural breaks method to find adequate group sizes for good visual interpretation. Afterwards I adjusted slightly the group-size to have nicer breaks.

So if for example breaks from jenks classification would have been:

1-12
13-21
22-42
43-50

than I could change them manually to:

1-10
11-20
21-40
41-50

by just clicking into the corresponding rows and changing the values.
Anything similar in Qgis?

Besides: Has anybody idea why Qgis always finds different group sizes any time I push on Jenks?


texts: how feasible is building my own map server? providing GeoCoding, Directions and Navigation

texts:
My need is a heavy use of routing. same as Google's Directions API. I'm looking for a viable alternative to abusing google and paying them too much.

how difficult is it to have my own server for that task?


texts: How to create a sector with different azimuths of one point using QGIS?

texts:
If I got one point like this:

ID longitude    latitude    azimuth

A   113.38107   22.35215    80

A   113.38107   22.35215    180

A   113.38107   22.35215    350

How to use QGIS to create a sector with different azimuths of one point like this?


texts: Using categorized and graduated visualizations with narrow table

texts:
I'm working with county data in a narrow table format with a WKT column. To simplify, if the headers are id, geo_level, metric, metric_unit, and metric_value, some entries for the county with id 5 are:

5 | county |  area  | square km | 500
5 | county | pop_married | percent | .7
5 | county | total_pop | people | 4,000
5 | county | per_immigrant | percent | .04
5 | county | households | residences | 2,900

Each county has twenty-eight rows, one for each of the twenty-eight metrics I'm tracking across all counties. So when I load it as a delimited text layer, QGIS renders 28 polygons for each county, stacked exactly on top of one another.

I'd like to be able to visualize specific metrics for each county ("here's % married from county to county", "here's population from county to county"). But visualizing with categorized or graduated symbols in the layer properties doesn't work, because the values in the metric_values are all in different units.

Is there a way to specify what row you'd like the visualization to come from? Or to specify based on a specific string in the desired row?


texts:
I am working with mapserver and c#
can anybody plz explain what does these line mean

function AddLayers()
{
var districtBase = new OpenLayers.Layer.MapServer(
        "District2",
        app.mapserverUrl() + "District2.map",
        { format: "OpenLayers.Format.Vector" },
        { isBaseLayer: true, singleTile: true, ratio: 1.5, transparent: 'true', opacity: .8}
    );
}

District2.map means a mapfile in mapserver?

app.mapserverUrl() + "District2.map means what?


texts: Landsat-8 and Sentinel-2: What are other options?

texts:
The data from Landsat-8 and Sentinel-2A&B satellite missions is available at no cost.
I am curious to know if there are other sources of freely available multi-spectral moderate resolution satellite imagery?
I am obliged for pertinent information about currently operational satellite missions.
EDIT:
The moderate resolution satellite imagery within spatial resolution of 10-30 meters, with regional to global coverage. I think imagery from Sentinel-2 and Landsat-8 falls in the category of moderate resolution.


texts:
Here is what I would like to do.  I have a polygon feature class called Subdivisions and a point feature class called Inspection.  If I label all the Subdivisions my map becomes too congested.  I only want to label Subdivisions that spatially contains an Inspection point.  If a Subdivision does not contain an Inspection point it does not get labeled.  Is there a simple way to do this?

The data is not setup to do joins or relates using primary keys. I know I can run a spatial join to accomplish this task but it involves extra steps an creating an extra feature class.  I have to maintain this map regularly so the few steps the better.

Is it possible to do this in ArcMap using python in label expressions or definition queries?  Let me know if you have a good solution.


texts: Preventing scaling during raster rendering in MapServer

texts:
I have encoded some metrics in a 4-band geotiff that will be served by MapServer (v7.4), and I want to "decode" them client-side for display using raster operations in OpenLayers (v6.1). (This is something somewhat similar to what is done in this openlayers workshop.)
The problem is that my metrics values range from 1 to 156 (and I used 255 for nodata), and MapServer scale them to 0-254 during rendering for better display, so the client do not have access to the original metrics.
I would need the client to have access to the original values. I have tried adding the following directives to the mapfile layer definition, but to no avail.
'PROCESSING "LUT=1:1,156:156"'
or 'PROCESSING "SCALE=0,254"'
Here is the layer definition:
LAYER
    NAME tif_source
    STATUS ON
    TYPE RASTER
    DATA "/test.tif"
    PROJECTION
        "init=epsg:32198"
    END

    METADATA
        wms_enable_request "*"
        wms_title "tif_source"
    END

    # PROCESSING "SCALE_1=0,254"  ## does not produce the expected result
    # PROCESSING "LUT=1:1,156:156" ## does not work either
    PROCESSING "LOAD_FULL_RES_IMAGE=YES"

    DEBUG 4
END

How should I prevent MapServer from scaling the raster values during rendering?


texts: Databinding ComboBox to GraphicsCollection - ESRI, WPF runtime

texts:
I'd like to bind a GraphicsCollection (which is a property in a view model) to a combo box in a XAML window. The problem is I can't figure out how to display attribute values in the combo box. I can get it to display the Geometry property like this:

<ComboBox DataContext="{StaticResource MyViewModel}"
                DisplayMemberPath="Geometry"
                ItemsSource="{Binding CollectionOfGraphics}"  />

But I want to be able to do something like this that will display the OBJECTID or some other attribute:

<ComboBox DataContext="{StaticResource MyViewModel}"
                DisplayMemberPath="Attributes[&quot;OBJECTID&quot;]"
                ItemsSource="{Binding CollectionOfGraphics}"  />

I also tried this approach - unsuccessfully:

<ComboBox DataContext="{StaticResource MyViewModel}"
                DisplayMemberPath="Attributes[0]"
                ItemsSource="{Binding CollectionOfGraphics}"  />

Platform: ArcGIS 10.1, WPF Runtime 10.1.1

Solution based on Antti Kajanus's answer below:

<ComboBox DataContext="{StaticResource MyViewModel}"
                ItemsSource="{Binding CollectionOfGraphics}" >
     <ComboBox.ItemTemplate>
          <DataTemplate>
              <TextBlock Text="{Binding Attributes[OBJECTID]}"/>
          </DataTemplate>
     </ComboBox.ItemTemplate>
</ComboBox>

Note that the attributite name (OBJECTID) did not have to be quoted.


texts: How to launch a QGIS core plugin from within a python plugin?

texts:
In this case I want to launch the Georeferencer core plugin from within a plugin I have written in python. But I don't understand how to interact with the core C++ plugin system from python, if that's possible. My plugin displays a QFileDialog and I optionally want to launch the Georeferencer plugin after looking at the selected files.

class ImportRasterDialog(QtGui.QFileDialog, FORM_CLASS):
    # ...

    def on_ok_clicked(self):

        files = self.selectedFiles()
        # Do stuff with the selected files
        # Somehow launch the Georeferencer plugin with a selected file

My searches keep coming back to this thread on the QGIS-developer mailing list, but this seems to be about lower-level access the functions of the Georeferencer plugin, whereas I only need to launch it (with a selected file): https://www.mail-archive.com/qgis-developer@lists.osgeo.org/msg22060.html

I couldn't see any suitable functions available from QgsInterface.

I am quite new to QGIS development and to this stackexchange, please excuse me if this question is obvious.


texts:
BACKSTORY:

I am currently in college (a nursing major - nothing to do with GIS or engineering), but I work for a local government's engineering department during the summers.  After a few changes with employees, I was assigned to keep the township's GIS mapping up to date.  I have been able to teach myself enough about ArcGIS Desktop to be able to create and update features in ArcMap, link to existing servers, and publish maps for ArcReader.

PROBLEM:

The local government is split into two different offices, with both needing the same maps.  We need a way for both offices to be able to update the maps without jeopardizing the other's work.  As of right now, we are only using one of our subscriptions, with employees from both offices using one computer to update their work and then using an external hard drive to transfer the files that are hyperlinked and the projects to the other office's server.  It is time consuming and a waste of our second subscription.

ATTEMPTED SOLUTIONS:

Since we have two subscriptions to ArcGIS Desktop (one for each office), we have decided that using our complimentary ArcGIS Online service is our best bet.  I have been able to create a web map that automatically opens to our area and connect to a server to add the parcel information, but I am unsure of the best way to connect our information from our desktop to online.  All of our layers are in four geodatabases. One solution may be to convert our geodatabases to a file, such as .csv, that is accepted by ArcGIS Online.  We have also attempted without success to connect our accounts using a portal or server, but we are not 100% sure we are doing this correctly.  I have been able to create a layer package in ArcMap and upload it to ArcGIS Online, but when I go to open the layer package, it only gives me an option to open it in ArcGIS Desktop, which doesn't help at all.

No employees that currently work for our organization have formal GIS training, though we do pay for a company to help us.  Unfortunately, the employee assigned to us will not offer any tips and only attempts to sell us a different online GIS system (but that system will not work for our needs/the cost would be far too much for the amount of information we would need to add to our online map).

How do we get our current ArcMap documents into ArcGIS Online?


texts: Reducing number of if statements in Field Calculator code block of ArcGIS Pro

texts:
I am using Python 3 in ArcGIS Pro. In plain English, my goal is calculate one field in a table using a list of values from another field. I can do this with a long list of if statements, but I can't figure out the syntax to do it with lists (and maybe for loops?).
For example, let's say I have two lists:
veg = ['Carrot', 'Celery', 'Broccoli']
fruit = ['Apple', 'Banana', 'Grape']

I would essentially like to ask..."If the field contains a value from the list 'veg', return 'Vegetable', else if the field contains a value from list 'fruit', return 'Fruit', else leave null"
I know this is possible I am just having trouble with the syntax. I would like to avoid doing this ad nauseum:
if FIELD == 'Carrot':
return 'Vegetable
elif FIELD == 'Celery':
return 'Vegetable'

I know it's possible I am just struggling with the syntax. I am reducing ~40 values in one field to 2 values in another.


texts: Getting query output from ogr2ogr when loading data?

texts:
I am trying to load a shapefile into a PostGIS instance using ogr2ogr. I am getting an odd error that I think might have to do with the specific version of PostGIS that we are using.

In order to debug what ogr2ogr is trying to do I would like to see the sql query that is generated behind the scenes.

I believe that Getting query output from ogr2ogr? is trying to get at the same thing. Like this OP, I am interested in reproducing the functionality available with shp2pgsql, which produces the create and insert statements as outputs with:

shp2pgsql -s 4326 shapefile.shp geodata_tablename > statements.sql

statements.sql then contains the create table statements that were used as well as the individual insert statements.


texts: Generating Map Tiles for MGRS Grid Graticule?

texts:
I'm possibly a bit out of my league here, but ready to get my hands dirty.  I want to generate map tiles for a MGRS grid graticule similar to http://tiles.mapbox.com/mapbox/map/afghanistan-mgrs.

I need to create a function that takes a map boundary (probably something like two lat,lng pairs for the map corners) and the output image dimensions and returns an image.  I'm looking at using Proj4 to convert map coordinates to pixel coordinates.  I'm looking for some open source or starter code for calculating the MGRS grid lines.  Most libraries I've seen are for translating a single coordinate which is fine for calculating the grid cell, but I'm not sure how to get the grid cell boundaries/plot the grid lines.

I'm looking at this for a mobile application so something in Java or C/Objective-C is my end goal.

Can someone give me a push in the right direction?

UPDATE

Refining my requirements a bit here...

I need to create a function that returns a 256x256 image given lat,lng boundary and zoom level.  The zoom is a standard slippy map like zoom where the lowest zoom level is 1 which as 1 256x256 tile for the entire world and zoom level 2 has 4 tiles, and so on.  For each tile and corresponding zoom level, I would like to draw the MGRS grid graticule boundaries and add a MGRS center label at the center of the MGRS grid cell.


texts: Automating/Improving workflow for adding HTML links to shapefile attribute tables in ArcMap

texts:
I'm doing some work that involves classifying aerial imagery. I'd like the consumers of this data to be able to access the photo associated with points through a hyperlink, to their location on a cloud storage service. I found good documentation on this process here:

http://resources.arcgis.com/en/help/main/10.2/index.html#//00s500000004000000

I've been able to add links one-by-one to shapefile attributes but there are many thousands of pictures in this dataset, and would take very long. I'm wondering if there might be some way to automate this process with a script of some sort, or perform a bulk job.
.


texts: Horizontal distance from point to polyline in ArcGIS Desktop?

texts:
I am trying to find a way to calculate the horizontal distance from a point to a polyline. I am not looking the shortest distance. I have a point layer with over 10,000 points and line layer. I am looking the horizontal distance. In another word. I am looking to calculate the distance from a point to polyline at the latitude of point in polyline.

You can see in the picture. I have a green line. I have also green points. I want to get those horizontal distances showed on red. I have some point on the left side of the line and some points on the right side of the line. I want to get horizontal distance (red line) to the line.

I tried to use near tool but it gives shortest distance. I tried use generate near table with angle and tried to find the distance with appropriate angle but, my polyline has a lot of curves so I am not getting the one in same latitude. I am using ArcGIS 10.3.1


texts: How to calculate/update multiple fields at once using QGIS 3.10 field calculator?

texts:
In QGIS 3.10, i want to calculate multiple fields at once on more than 10 similar layers.
I now i can do it manually, but I'm looking for an automatic way.
I need something like this:
"nomrue"  should be filled with "potable_nomrue" And

"debut"  should be filled with "potable_debut" And

"fin"  should be filled with "potable_fin" And

"seqsegment"  should be filled with "potable_seqsegment" And

"idnoeuddebut"  should be filled with "potable_idnoeuddebut" And

"idnoeudfin"  should be filled with "potable_idnoeudfin" And

"vocation"  should be filled with "potable_vocation" And

And so on...
And i need to repeat the operation on multiple layers.
I found some Python code, but i don't understand how to use Python...
Ihave:
lyr = iface.activeLayer()

with edit(lyr):
    for f in lyr.getFeatures():
            f['nomrue'] = f['potable_nomrue']
            f['debut'] = f['potable_debut']
            f['fin'] = f['potable_fin']
            f['seqsegment'] = f['potable_seqsegment']
            f['idnoeuddebut'] = f['potable_idnoeuddebut']
            f['idnoeudfin'] = f['potable_idnoeudfin']
            f['vocation'] = f['potable_vocation']
            f['anneerehab'] = f['potable_anneerehab']
        lyr.updateFeature(f)

But there's an error... What's wrong?


texts: ArcMap lagging when I use the measure tool

texts: ArcMap is lagging when I click the measure tool, my data isn't particularly big (biggest .shp is ~150mb) and my PC is brand new. Any idea what may be causing this?


texts: Arcpy cannot locate rasters inside geodatabase

texts:
I've bunch of rasters inside a geodatabase and I'm trying to call them as inputs in a script, but I'm running into a strange problem. The script cannot find/locate these rasters, so I manually added and inspected the rasters. They appear to be valid and not corrupt. Next, I exported these rasters as TIF files and tried locating these files. Here are the results:

# Trying to call in rasters from my gdb
import arcpy, os

coName = "KENT_10001"

DEV = os.path.join("C:/Imagery/", str(coName) + "_Masks.gdb/", str(coName) + "_DEV")

os.path.exists(DEV)
False
os.path.isfile(DEV)
False

#Trying to call in rasters exported from gdb as TIFs
import arcpy, os

coName = "KENT_10001"

DEV_TIF = os.path.join("C:/Imagery/", str(coName) + "_DEV.tif")

os.path.isfile(DEV_TIF)
True
os.path.exists(DEV_TIF)
True

I don't believe the geodatabase is corrupt, because we were successful in running other models using the same inputs from the geodatabase. Additionally, I've tried executing the same test using explicit paths as well instead of using os.path.join and the results are still false for both path and file.

Update

arcpy.Exists(DEV) returns True value, but when I try to call this raster into raster calculator. It gives me the following error:

File "<string>", line 5, in rcexec NameError: name 'DEV' is not defined Failed to execute (RasterCalculator)


texts:
I recently found an issue with Geotools and getting the datastore for a shapefile.

If I try:

FileDataStore dataStore = FileDataStoreFinder.getDataStore(file);

-or-

Map connect = new HashMap();
connect.put("url", file.toURL());

DataStore dataStore = DataStoreFinder.getDataStore(connect);

getDataStore() always returns null. But using:

ShapefileDataStore dataStore = new ShapefileDataStore(file.toURI().toURL());

then the dataStore is set to the shapefile. Anyone have any ideas why this is? Is this a bug with Geootools 8.2? I'm running Geotools on Tomcat 7, perhaps there is some issue with Tomcat and the DataStore and FileDataStore classes?

Thanks!


texts:
I am using ArcMap 10.4 and I want to convert my layer of Athens (http://geodata.gov.gr/dataset/api0unon-o1kodou1kwv-tetpaywvwv-anuou-a0nvaiwv/resource/d2e17dd6-8e1e-4c34-9111-fb3ac1bc9da4) municipality which is in

Current coordinate system:

GGRS87_Greek_Grid
Authority: Custom

Projection: Transverse_Mercator
false_easting: 500000,0
false_northing: 0,0
central_meridian: 24,0
scale_factor: 0,9996
latitude_of_origin: 0,0
Linear Unit: Meter (1,0)

To

WGS_1984_Web_Mercator_Auxiliary_Sphere
WKID: 3857 Authority: EPSG

Projection: Mercator_Auxiliary_Sphere
False_Easting: 0,0
False_Northing: 0,0
Central_Meridian: 0,0
Standard_Parallel_1: 0,0
Auxiliary_Sphere_Type: 0,0
Linear Unit: Meter (1,0)

My goal is to make a web map and use also OpenStreetMap as a basemap.

I have tried to use the project tool but it does not have this transformation.

I have also tried to make a new transformation as described in this link https://georepository.com/transformation_1272/GGRS87-to-WGS-84-1.html

Parameter   Value / File Ref    Sign Reversible Unit
X-axis translation  -199.87 YES metreOpen
Y-axis translation  74.79   YES metreOpen
Z-axis translation  246.62  YES metreOpen

and still I get an error.


texts: PostGIS polygon edge analysis (orientation, edge length)

texts:
I am rather new to the world of GIS and especially PostGIS, so please excuse me if the answer seems evident...

I would like to do analysis on a number of buildings. One thing I am interested in is their facade surfaces along with the respective orientation. As illustrated in the picture below, I would like to have the length and (normal) orientation of all edges in a series of polygons. In the example I highlighted only one surface.

A result table could look like this:

building_id | edge_id | orientation | edge_length
-------------------------------------------------
      1     |    1    |     315     |    10.0
      1     |    2    |      45     |     7.0
      1     |   ...   |     ...     |     ...

However, I am not sure if it is a smart way to store the result for further processing (e.g. calculate distance from edge to next building, etc.). So my question is twofold:

Is there an efficient PostGIS function that can analyze a polygon's edges? In case there is no native PostGIS function I would alternatively be interested in an Python-based approach.
What would be a smart way to store the result in a PostGIS table, since the polygons may have different numbers of edges?


texts: Loading simple ESRI HTML JavaScript in Internet Explorer 11?

texts:
Is there any particular reason this simple ESRI HTML JavaScript won't load in IE11? Chrome works fine.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="http://js.arcgis.com/3.12/esri/css/esri.css">


texts: How do I edit the rules in an already created geodatabase topology in ArcGIS (10.7.1)?

texts:
I'm using ArcGIS Desktop 10.7.1 with an Advanced license and created a topology for a file geodatabase I'm working with. (I have complete access to the data.) I'd like to edit the rules that I've created, but when I try to edit them via the Rules tab in the Layer Properties dialog for the topology layer in ArcMap, all of the options (Add Rule, Remove Rule, etc.) are greyed out except for Save Rules and Description. (See attached screenshot).

Is there a way to edit the geodatabase topology without removing it and starting a new one?


texts: OpenLayers 3: Cross-Origin Request Blocked: The Same Origin Policy disallows

texts:
Using OpenLayers 3, I cannot get this message to go away:

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON. This can be fixed by moving the resource to the same domain or enabling CORS.

This is the code:

// Ol3 only supports Projections "EPSG:4326" and "EPSG:3857". For every other projection you need proj4js
        proj4.defs("EPSG:2236", "+proj=tmerc +lat_0=24.33333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs");

        // Leases Layer
        var myLayer = new ol.layer.Vector({
            source: new ol.source.GeoJSON({
                projection: 'EPSG:2236',
                url: 'http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON',
                crossOrigin: null
            })
        });

// View
        var view = new ol.View({
            projection: 'EPSG:2236',
            center: [0, 0],
            zoom: 4
        });

        // Map
        var map = new ol.Map({
            target: 'map',
            renderer: 'canvas',
            layers: [myLayer],
            view: view
        });

I have tried setting the crossOrigin setting to:

crossOrigin: null
crossOrigin: 'null'
crossOrigin: 'anonymous'

I only see the zoom in/out control but the layer is not rendered.

I went with simon's option 3 below.  I enabled CORS in GeoServer by copying the necessary jetty-servlets jar files and enabling it in the \WEB-INF\web.xml:

<filter>
    <filter-name>cross-origin</filter-name>
    <filter-class>org.eclipse.jetty.servlets.CrossOriginFilter</filter-class>
    <init-param>
        <param-name>allowedOrigins</param-name>
        <param-value>*</param-value>
    </init-param>
    <init-param>
        <param-name>allowedMethods</param-name>
        <param-value>*</param-value>
    </init-param>
    <init-param>
        <param-name>allowedHeaders</param-name>
        <param-value>*</param-value>
    </init-param>
</filter>
<filter-mapping>
    <filter-name>cross-origin</filter-name>
    <filter-pattern>/*</filter-pattern>
</filter-mapping>

After I did that, I tested the page again and receive the same error:

Cross-Origin Request Blocked: The Same Origin Policy disallows reading the remote resource at http://myserver:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view1&maxFeatures=1&outputFormat=JSON. This can be fixed by moving the resource to the same domain or enabling CORS.

Looks like I am still missing something.  Do I have to do anything from the OpenLayers Side?

I ended up getting rid of Jetty and uninstalling GeoServer completely.  The problem is when you install the geoserver windows installer, it installs a version of jetty that is 4 years old! (Jetty version 6.1.8)  Even though I had copied the jar files for CORS, it is only supported in Jetty 7+.

I found out that you can install a WAR file.  I decided to use Tomcat since that is what GeoServer is mostly tested on according to this note from the GeoServer website:

Note GeoServer has been mostly tested using Tomcat, and therefore these instructions may not work with other container applications.

These are the instructions for installing the WAR file:

http://docs.geoserver.org/stable/en/user/installation/war.html

This is a nice how-to video also:

https://www.youtube.com/watch?v=YEOA8WWWVCw

After you complete the install, you then enable CORS:

http://enable-cors.org/server_tomcat.html


texts: How to import for QgsVectorLayerJoinInfo in QGIS3 (python)?

texts:
I am trying to join layers programmatically following the answers in the GIS.SE question -
Join table field with shapefile programatically (via PyQgis)? :

# Get input (csv) and target (Shapefile) layers
shp=iface.activeLayer()
csv=iface.mapCanvas().layers()[0]

# Set properties for the join
shpField='code'
csvField='codigo'
joinObject = QgsVectorJoinInfo()
#...

But I get an error about QgsVectorJoinInfo not defined. I tried (unsuccessfully) to import it with:

from qgis.core import QgsVectorJoinInfo

(because the documentation linked above says that QgsVectorJoinInfo is in the "core" library).

How should I import it properly or otherwise make the above code work?

(In general, how to tell from documentation which library to import?)


texts: Why is intersect giving blank layer in QGIS?

texts:
I have two shapefiles. One is broken into the regions I need to use, the other has much more detail along the coast line. I'd like to use either intersect or clip to get a more detailed coastline on the eastern edge of the map, but I wind up with a blank layer with both.

With vector > geoprocessing tools > intersect the new layer is empty no matter how I arrange the settings. If I use vector > geoprocessing tools > clip and put the big unified polygon as the "input vector" I wind up with the opposite of what I want -- a new layer that has one big polygon and no coastal detail. What I want is a layer of small polygons with coastal detail. What am I doing wrong here?

I'm working in QGIS, though if this would be easier to do in PostGIS I can also use CartoDB.


texts: Simplification Tolerance value field is missing from Simplify Polygon tool of ArcGIS for Desktop?

texts:
Trying to run Simplify Polygon (Cartography -> Generalization) tool in ArcGIS 10.2 (ArcInfo License).

I ran the tool successfully once, then tried again on a new polygon, but the simplification tolerance value field is missing...so running the tool results in an error that the missing value is required.

I tried restarting ArcMap but the field is still missing.

Is this a bug in the software or has anyone else had a similar experience?


texts: How can ArcGIS Editor for OpenStreetMap preserve Unicode data?

texts:
Currently i am using ArcGIS Editor for open Street map for conversion of OSM file to file Geodatabase. But after conversion geometry is perfect but there is no Unicode data available in the attribute table. Am i missing some thing?

I am particularly interested in the house numbers and buildings names.


texts: Can't connect to foreign SQL Server table from PosgreSQL using OGR_FDW

texts:
I'm using Windows 7, with PostgreSQL 9.3.6 installed and PostGIS 2.1.5. While I wish I could get everything into PostgreSQL, I can't and need to connect to another server running SQL Server. I need to use a foreign data wrapper for this, and I've installed ogr_fdw as described here.
I'm trying to connect a non-spatial table in the foreign SQL Server to my PostgreSQL, and will then join that foreign table to a PostgreSQL spatial table.

From the instructions, I'm unclear how to setup the parameters correctly to make the connection (assuming it is even possible to get it working), and I don't understand how to use the ogr_fdw_info utility. For example, I tried:

ogr_fdw_info -s dbname='databasename' host='servername' user='SQL-serverusername'

but that doesn't get me anywhere. Does anyone know how to connect to a foreign SQL server table using this ogr_fdw library?

With Paul's suggestion, I tried using ogrinfo like this:

ogrinfo "MSSQL:server=ServerName;database=DatabaseName;UID=SQLServerUserName;PWD=SQLServerUserNamePassword"

When I tried that, at least I got

Error initializing the metadata tables: [Microsoft][ODBC Server Driver][SQL Server]CREATE TABLE permission denied in database 'DatabaseName'.

To me that means at least I'm sort of connecting, but not completely. The SQL Server Admin gave me a username with only read rights, and that seems logical since I only want to read one non-spatial table and put it as a foreign table inside my PostgreSQL database on my own machine. I'll have to make sure I can connect using ogrinfo before I can begin to understand how to use ogr_fdw. I tried adding "estimatedmetadata=true" at the end of my ogrinfo statement, but that didn't make a difference. I don't understand exactly why the error would be in creating a table, since I'm only wanting to connect, not create. Any additional hints on getting this connection to work?

With Paul's suggestion, I added the environment variable to my Windows machine:

MSSQLSPATIAL_USE_GEOMETRY_COLUMNS=No

and now the error has gone away. Then I also found I needed to add a second environment variable:

MSSQLSPATIAL_LIST_ALL_TABLES=YES

This enabled me to correctly connect to the SQL server using OGRINFO and see the non-spatial tables. Now however, when I try:

ogr_fdw_info -s "MSSQL:server=ServerName;database=DatabaseName;UID=SQLServerUserName;PWD=SQLServerUserNamePassword"

It says it can't connect to that data source. I'm puzzled since the connection string works fine with OGRINFO, but not OGR_FDW_INFO.


texts: Running QGIS processing tools in the OSGEO4W shell using the qgis_process functionality, syntax errors

texts:
I am trying to run QGIS processing tools in the OSGEO4W shell using the qgis_process functionality. I believe this functionality is fairly new as I have been having a hard time finding documentation for the syntax.
I am trying to run a raster calculator command using two tiffs as input, everything about the command runs fine except for the LAYERS parameter always returns an error as I am not sure how to list two layers in the parameter.
the command is:
qgis_process run qgis:rastercalculator --CELLSIZE=0 --CRS="ESPG:32610" --EXPRESSION=""octclip@1" + "marclip@1"" --EXTENT="D:/UPCD/clip1/octclip.tif" --LAYERS=D:/UPCD/clip1/octclip.tif D:/UPCD/clip1/marclip.tif --OUTPUT="D:/UPCD/clip1/try1.tif"

I have tried many things such as
--LAYERS=D:/UPCD/clip1/octclip.tif D:/UPCD/clip1/marclip.tif

--LAYERS="D:/UPCD/clip1/octclip.tif" "D:/UPCD/clip1/marclip.tif"

--LAYERS="D:/UPCD/clip1/octclip.tif","D:/UPCD/clip1/marclip.tif"

--LAYERS="D:/UPCD/clip1/octclip.tif":"D:/UPCD/clip1/marclip.tif"

to no avail.
The documentation says it will allow paths both in quotes as well as without but I can not find anything regarding passing two paths to the parameter argument.
First time question poster, let me know if there is any information I am missing.


texts:
I am building a stand alone .NET applications using ArcObjects.
I want to run it on the server, where I have ArcGIS for Server.

I am trying to bind the license at runtime to the ArcGIS server runtime.
But the application giving an error message "Invalid ArcGIS Runtime Binding"

I am using this command

RuntimeManager.Bind(ProductCode.Server)

I don't want to install the engine there or have a license for the engine, because the application run at night when there is no activity , so I want to use the server license

Is that possible?


texts: Split line at a point with ArcGIS 10.1 Basic level license?

texts: I have a point and polyline shapefile and want to split lines using points -- want to emulate the precise functionality of Split Line at Point tool (Data Management). I have ArcGIS 10.1 Basic level license and don't want to use third party tools -- want to achieve this by using Basic level license geoprocessing tools available within ArcGIS 10.1. Any suggestions would be very much appreciative.


texts:
Yet another intersection issue. I'm working with OSM database in PostGIS.
Now trying to count all buildings that fall into a countries bounding box.

This is my query (bbox for Germany).

SELECT
  COUNT(osm_id) AS num
FROM
  planet_osm_polygon
WHERE
  tags->'building' IS NOT NULL AND
  way && ST_Transform(ST_MakeEnvelope(12.8719,52.3700,13.9382,52.6646, 4326), 3857)

No matter whether I filter for buildings or not, the query is slow.
My expectation would be below 10s but reality is about 100 seconds.

This is my query plan:

Aggregate  (cost=500537.88..500537.89 rows=1 width=8)
 -> Bitmap Heap Scan on planet_osm_polygon  (cost=7691.95..500254.94 rows=113173 width=8)
    Recheck Cond: (way && '0103000020110F0000010000000500000076B9815A3DDD354137FCF32B5A325A4176B9815A3DDD3541133B4058FB665A410E689A53E9AC3741133B4058FB665A410E689A53E9AC374137FCF32B5A325A4176B9815A3DDD354137FCF32B5A325A41'::geometry)
    Filter: ((tags -> 'building'::text) IS NOT NULL)
     -> Bitmap Index Scan on planet_osm_polygon_index  (cost=0.00..7663.66 rows=133748 width=0)
        Index Cond: (way && '0103000020110F0000010000000500000076B9815A3DDD354137FCF32B5A325A4176B9815A3DDD3541133B4058FB665A410E689A53E9AC3741133B4058FB665A410E689A53E9AC374137FCF32B5A325A4176B9815A3DDD354137FCF32B5A325A41'::geometry)

Eventually this should run against another table with 200+ country polygons.


texts: I want to do the following: I have a set of cells with the value 1, and a set of values with 2. I want to calculate for each cell [2] the distance - ideally number of cells - to the nearest cell 1. Can I somehow do this in GRASS GIS?


texts: Moving points closer to polyline using ModelBuilder?

texts:
I have a lot of points I need to move that are on the edge of the parcels closer to the centerline(road), not directly on top but about 5 feet away from the centerline.

I know i can edit and move the points manually but is there a way to automate the process using ModelBuilder?


texts: Getting attribute info from OpenLayers (QGIS2Web)

texts:
The below files are generated from QGIS2Web tool and altered by me, I generated the OpenLayers code using QGIS2Web. The problem is not showing the pop-up(attribute information) for multiple Layer files. For a single file pop-up function is working fine (I attached the screenshot with pop-up).

I would like to show my attribute details on the right side of the panel(permanent) when the user clicked the particular polygon.

How to customize based on my need?

The above pic shows the pop-up function (working fine for single layer).

The above pic shows the multiple layer structure. I want to show attribute information for the topmost visible layer.


texts:
I have set up  a Raspberry PI with a RasPiGNSS GPS-Module. After some binary convertion I get the following RINEX(ftp://igs.org/pub/data/format/rinex210.txt) files.

.nav
.obs
.sbs

I'm a little confused about the difference between the obs and sbs files.

What information contain these files and from where the gps-receiver gets the messages?

If my understanding is correct: I get the sbs from a base station (in my case i will be a EGNOS in Zürich) and it contains real time information about time and orbit corrections.


texts: Invalid Primary Key error message when attempting to join table view and feature layer

texts:
I am attempting to join a feature layer to a table view in order to export a table of all non-joining features in the table view by selecting NULLS from the join field related to the feature layer. When I attempt using AddJoin to accomplish this I am met with the following error.

  ERROR: 999999: Error executing function. The layer has an invalid
  primary key

I have used AddJoin to join these two datasets earlier in my script successfully, but that was joining the Table View to the Feature Layer, not the other way around. Please see below for more details and code.

I create the Table View and confirm it exists with this bit of code:

for file in gdb_files:
    print("EDes table available: {}".format(file))
    if 'GIS' in file:
        print("Creating Table View for EDes file")
        arcpy.MakeTableView_management(os.path.join(gdb_path, file), "EDES_TableView")
        print(arcpy.Exists("EDES_TableView"))

I create the Feature Layers and confirm they exists with this bit of code:

arcpy.MakeFeatureLayer_management(os.path.join(gdb_path, "Some_SDE_FC"), "Some_F_Layer")
arcpy.MakeFeatureLayer_management(os.path.join(gdb_path, "Some_SDE_FC_2"), "Some_F_Layer_2")

print(arcpy.Exists("Some_F_Layer"))
print(arcpy.Exists("Some_F_Layer_2"))

I am able to join a table view and feature layer successfully prior to the error message with this bit of code:

arcpy.env.workspace = gdb_path
arcpy.AddJoin_management("Some_F_Layer", "BBL", "EDES_TableView", "BBL", "KEEP_ALL")
arcpy.CopyFeatures_management("Some_F_Layer", "Match")

But when I attempt to rejoin the same Table View and Feature Layer, but in the opposite order, I get the above mentioned error. See code bit below:

arcpy.AddJoin_management("EDES_TableView", "BBL", "Some_F_Layer", "BBL", "KEEP_ALL")


texts: Cannot access `properties` key/value content when querying ArcGIS Feature Service for GeoJSON

texts:
Using ArcGIS with a developer account, I created a new layer, uploaded a .geojson file, and happily followed along the "query a feature layer" tutorial on their labs site. I successfully generated an OAuth Token in Postman, but when I changed the "f" (format) value in the x-www-form-urlencoded body to "geojson" it looks similar to the originally GeoJSON file, but the properties are missing/null and the only key info included are the Polygon coordinates. The data remains in the ArcGIS Dashboard but isn't included in the data. Any advice? Is this a bug?
Original File (.geojson):
{
  "type": "FeatureCollection",
  "features": [{
    "geometry": {
      "coordinates": [
        [
          [50.0115001901903, 50.2259915257836, 0],
          [50.0115418500254, 50.225991583611, 0],
          [50.0115447598285, 50.2242567612666, 0],
          [50.0115031005844, 50.2242567034437, 0],
          [50.0115001901903, 50.2259915257836, 0]
        ]
      ],
      "type": "Polygon"
    },
    "ID": "asdf98d0-fa90-4293-a780-a6988f5486f7",
    "id": 1,
    "type": "Feature",
    "properties": {
      "CreationDate": "Wed, 17 June 2020 18:40:41 GMT",
      "Editor": "dkoch"
    }
  }]
}

Received file from POST /0/query rest endpoint with f: geojson, token: (working token), where: 1=1:
{
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [
                            50.0115001901903,
                            50.2259915257836
                        ],
                        [
                            50.0115418500254,
                            50.225991583611
                        ],
                        [
                            50.0115447598285,
                            50.2242567612666
                        ],
                        [
                            50.0115031005844,
                            50.2242567034437
                        ],
                        [
                            50.0115001901903,
                            50.2259915257836
                        ]
                    ]
                ]
            },
            "properties": null
        }
    ]
}


texts:
I was given a point shapefile that does not have coordinate system. So, I need to specify it. I was told that this have following reference system.

Coord. Sys. – 17 South (UTM)
Units – US Survey ft
Datum – WGS 1984 (7P)
Geoid model – GEOID12A (Conus)

I have never seen options for specifying Geoid model in ArcGIS, and also, I have not seen WGS 1984 (7P). Except this, I am familiar with other information. These two things are confusing me. Will I be fine if I just use WGS_1984_UTM_Zone_17S?

Can someone help me if this is fine. Please suggest.


texts:
I have a raster and a vector point file.
I need to generate another raster of same resolution showing how many points are inside the corresponding pixels of the original raster.

Something like that:

Is there a way of doing that in QGIS?

I just want the number of points, no attributes involved.


texts: Why are short integer fields on SDE auto-updating to Precision 5?

texts:
I am having a perceived small issue in Server for ArcGIS (10.3.1 and 10.4.1 on SQL Server) where whenever a Short Integer field is added to a feature class, the Precision auto-updates to Precision: 5 regardless of configuration when created.

In the screenshots below, I have added a short integer field; InstallYear. InstallYear has a Precision of 4.

Once I save the changes and navigate back to the properties of the feature class, you can see that it has auto-updated upon save to Precision: 5. The field was added in ArcCatalog. I have also tested this by adding a field in ArcMap with the same result.

File geodatabases short integer precision cannot be altered. However, I tested it on a shapefile and the configured precision does save.

Can anyone explain why short integer precision is auto-updated to 5 on SDE, regardless of configuration? Referencing an Esri document, short integer precision of 1-5 is supposedly supported on SQL Server unless I am reading that incorrectly.


texts:
I'm a bit lost with the zonal statistics tool in QGIS. I have a raster layer representing population densities per km2 for Asia. I've clipped the raster to my Vietnam provinces polygon as it is my area of interest. For the analysis, I just need to know how much population is present inside each province. Now I've run the zonal statistics tool, and it gives results as presented below:

I just don't understand if the "popsum" field has taken all the pixels and corresponding values and is presenting this sum just for 1km2, or has it actually sum up the entire pixels for the entire area of the polygon? I've set my CRS to UTM Zone 48N corresponding to Vietnam and wider south east asia.


texts: I'm trying to rotate multiple features  in a certain angle, but each on its own pivot, is there a way to do it in ArcMap? Whenever I try to select the features and use the rotate tool all of them rotate on the same axis.


texts: ArcGIS Javascript Viewer with ArcGIS Server

texts:
I'm certain this question must have been asked/answered, but I just can't find it.

I'm running ArcGIS Server and want to build a basic viewer that points to my locally hosted services.  I'm not interested in doing anything with ArcGIS.com; however all the examples hosted at esri illustrate interfacing with ArcGIS.com.

How can I point to my local map?

Alternatively, can you point me to where I can find the answer (as I'm certain this must have been asked already)


texts:
SpatiaLite supports a remarkable number of different dimension models.

Interestingly, some of their values support a "measure" coordinate (XYM, XYZM). How would this measure-coordinate be used, and what is the correct term for searching to figure out how it is used?


texts:
Anyone know if there's some secret magic parameters that need to be set for the Extend and Trim (Editing Toolbox) tools in ArcMap to run properly?

I have a dataset of about 2000 "interior" lines which have been systematically separated from their surrounding "exterior" lines, and need to be extended or trimmed in order to reconnect with their borders.

I can do this manually one each line using Extend/Trim in the Advanced Editing toolbar, but as I need to do 2000+ records (plus a possible 16 000 for another project) I really need to do it in batch mode. This would mean using the Extend/ Trim tools from the Editing Toolbox. However, I can't seem to get these tools to actually run properly.

When I run Extend (Editing Toolbox) the tool runs 'successfully' in about 0.45 s with no error messages, but doesn't actually do anything. I have several test areas where I have measured the distance of the gap between the line that needs to be extended and the next line, these gaps are 0.5-0.8m maximum. I have set the Extend distance to various lengths such as 1.0m, 1.5m, 2.0m, as well as leaving it null, which should extend all lines to the next nearest feature. It still does nothing, the lines in my test areas do not appear to be altered at all. I've read through the ArcGIS documentation and a variety of forum posts, and I still can't figure out what I'm missing. I've had a similar lack of results when using the Trim tool.

Things I have already tried:
- running tool in and out of an Edit session
- running tool with lines to be extended selected or un-selected
- extending lines from one feature class to another
- extending lines from one feature class to itself
- running tool manually from toolbox
- running tool programmatically through Python window
- running FindIdentical (Editing Toolbox tools) and DeleteIdentical (Editing Toolbox) in case there are any overlapping features which are causing confusion
- running tool with and without other feature classes present in the view
- checking the geoprocessing and tool Environments to ensure that processing extent, projection, linear measurement units, etc. are logical

Anyone have any idea what I might be missing?

I am using ArcGIS for Desktop 10.2.2 with an Advanced license.

Here is an example of an area of "interior" (blue) lines that need to be extended to the "exterior" (green) lines.


texts: Deleting feature class features faster using ArcPy?

texts:
I have written this script that recurses through all feature classes in a given set of feature datasets and deletes all their features. It appears to do the job but it runs pretty slowly. Is there anything i'm doing wrong here or are there any obvious ways to speed things up?

I'm using DeleteFeatures_management to do the deed. DeleteRows_management also seems to work.

import sys
import os
import arcpy
from arcpy import env
import datetime
import getpass

try:
    passwd = getpass.getpass("Enter the sde user password: ")

    sdeConnectionFileDir = os.environ.get("TEMP")
    databaseName = ""
    fileName = "temp.sde"

    # Delete any pre-existing SDE connection file.
    fullPath = sdeConnectionFileDir + '\\' + fileName
    if os.path.exists(fullPath):
        os.remove(fullPath)

    # Create temporary SDE connection file.
    arcpy.CreateArcSDEConnectionFile_management (
        sdeConnectionFileDir, fileName,
        "sdeserver", "5151", "",
        "DATABASE_AUTH", "my_sde_user", passwd,
        "SAVE_USERNAME", "SDE.DEFAULT", "SAVE_VERSION"
    )

    env.workspace = fullPath

    # ArcPy status codes.
    returnCodes = {'WARN' : 0, 'INFO' : 1, 'ERROR' : 2}

    featureDatasets = []
    featureDatasets.extend(arcpy.ListDatasets("dataset1*"))
    featureDatasets.extend(arcpy.ListDatasets("dataset2*"))
    featureDatasets.extend(arcpy.ListDatasets("dataset3*"))

    list = '[%s]' % ', '.join(map(str, featureDatasets))
    response = raw_input("\n***** WARNING!!! ***** \nAll data will be deleted from all feature classess in the following datasets: \n\n" + list + "\n\n |--> Type DELETE to begin removal: ")
    if response == "DELETE":
        print "\nStarted: " + str(datetime.datetime.now()) + "\n"
        for dataset in featureDatasets:
            print "Processing dataset: " + dataset
            for fc in arcpy.ListFeatureClasses("*", "ALL", dataset):
                rowCount = int(arcpy.GetCount_management(fc).getOutput(0))

                if rowCount > 0:
                    print "  -- Processing feature class: " + str(fc) + " (" + str(rowCount) + " rows)"
                    #arcpy.DeleteRows_management(fc)
                    arcpy.DeleteFeatures_management(fc)

        print "\nCompleted: " + str(datetime.datetime.now())

except Exception as e:
    if arcpy:
        arcpyErrors = arcpy.gp.getMessages(returnCodes['ERROR'])
        if arcpyErrors:
            sys.stderr.write(arcpyErrors + "\n")
    sys.stderr.write(str(e) + "\n")
    sys.exit(1)

EDIT --

I put some performance timers in the script and here is the data:

Time to retrieve datasets: 0:00:01.254000
Total Feature Classes: 1682
Total Feature Classes with Data: 124
Total Features Processed: 190222
Total Run Time: 3 hours, 16 minutes

The breakdown:

Feature dataset --> list feature class calls:

* AVG   0:00:02
* MIN   0:00:01
* MAX   0:00:07
* COUNT 40
* TOTAL 0:01:08

Feature count calls (majority of the time):

* AVG   0:00:06
* MIN   0:00:01
* MAX   0:00:16
* COUNT 1682
* TOTAL 2:41:00

Feature deletion calls (reduced because only feature classes with rows are processed):

* AVG   0:00:17
* MIN   0:00:02
* MAX   0:03:22
* COUNT 124
* TOTAL 0:34:31


texts:
I am wondering if anyone knows if it is possible to extract the polygons shown in the following embedded Google Map. I cannot seem to find anything using 'Inspect' or 'View Page Source'
https://sustainability.aprilasia.com/en/


texts: How to clean layers from non topological databases?

texts:
I'm doing my internship in a consulting office which is switching from MapInfo to QGIS/GRASS 1.8.0.

I posted here a few days ago about the way to clean a polyline from a Lidar data with v.clean.

I've tried all algorithms possible, made every test I could try, without obtaining a cleaned polyline without any topological error.

Now I'm working on many other projects, and I realize that most of data coming from MapInfo or ArcGIS users are full of topological errors for QGIS/GRASS.

I'm feeling fear growing up, because I still didn't find the way to clean all those data with QGIS/GRASS. And I've tried all I could.

I used break then rmdangle to erase dangles with a good threshold, but it didn't work every time. I tried rmarea to erase some little islands of lines but it don't work if you want to keep the main polyline, I tried snap to full the holes in the line but it creates more artifacts if the threshold is too high, and it don't snap if the threshold is too small.

I think that it is obviously a big question that anyone working with QGIS/GRASS in a professional environment must face. But on the other hand I'm very surprised to see that many people like me asking help for this problem finally never find answers. Is it a taboo subject ? Or is it so simple to use v.clean that there's no comment to do ?

So I'm addressing those who often need at work to clean their layers with a hundred errors of all kind (double geometry, dangles, self-intersection...).

How do you success at cleaning it until you have 0 errors ?

Usually, people says "v.clean is powerful, you just need to launch rmdangle for dangles...". I just give you a quote from GRASS documentation : " Topological errors can be corrected either manually within v.digit or, to some extent, automatically in v.clean."

What is "some extent" ? From my small experience, it seems not to be a very large extent...and if this extent is very small, do you correct all the other errors manually ? Is there more option than v.clean to correct topology ?

I promise that if I find the answer, I'll post a very complete tutorial for beginners and idiots like me.

Thank you very much for your help, I think we are numerous to struggle with topology !!


texts: ArcObjects: Setting image of a command button to a partially transparent bitmap

texts:
In an ArcMap 10 add-in, I am attempting to create a toggle button of sorts whose image changes based on the button's state.

This answer describes how to set an add-in button's image programmatically and I am using the same technique. The problem is that this method does not appear to support full alpha channel transparency, only single-color transparency, so those fancy PNG icons they added at ArcGIS 10 don't look quite right.

The button's image is set initially by the add-in framework (using a 32-bit PNG file with alpha channel) -- the image looks as you would expect with proper partial transparency.

To create the image toggling behavior, in the button's OnClick method, I check the state of a boolean variable and set the button's image (using the technique in the above-mentioned answer) to one of two bitmaps stored as resources in the add-in assembly. This all works fine except that the button's image is given a gray background instead of a transparent background.

I have searched around a bit and come to the conclusion that ESRI's GetIPictureDispFromBitmap method converts from a GDI+ bitmap to a GDI bitmap and in doing so loses the alpha channel of the original bitmap.

I have found a solution to this issue and will post it as an answer shortly.


texts: Georeferencing a raster on top of 2 or more polygons

texts: I am trying to georeference raster sketch maps that have polygons drawn on them to line up with vector polygons.  Its relatively easy to do this using the georeferncing toolbar in arcmap for 1 polygon, however once there are 2 polygons it becomes much more difficult to get it accurate.  Does anyone have any experience with this?  See image attached.


texts: Method to split line segments in a road network that are greater than X length into n smaller line segments using ArcGIS?

texts:
I'm using ArcGIS Desktop 10.2 and I'd like to take roads in a network greater than a certain length and split them into equal smaller pieces. I've noticed that splitting by percent can be achieved via the Editor Menu; however, it appears this must be done on a line-by-line basis.

Could you suggest ways of automating the process to split all roads over a certain length into smaller parts based on the percent of the original?


texts: Band unit change from mol/m^2 to micromol/m^2 in GEE

texts:
I am working on NO2 data of Sentinel 5p in GEE(Google Earth Engine) the unit of band "NO2_column_number_density" is  mol/m^2.
Will someone guide me how can I change the unit from  mol/m^2 to micro mol/m^2?
var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
  .select('NO2_column_number_density')
  .filterDate('2019-06-01', '2019-06-06');

var band_viz = {
  min: 0,
  max: 0.0002,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};

Map.addLayer(collection.mean(), band_viz, 'S5P N02');
Map.setCenter(65.27, 24.11, 4);


texts: How realistic is it to be able to build a dynamic map with ~50k points that change colour depending on supplied data?

texts:
I'm a programmer, but I'm new to GIS.

I want to build a map that displays the status of some IoT devices I work with.

My ideal vision is a world map, which would fill the screen, but for now, Europe would do, or even a single European country. It would be nice to be able to switch between EU and a single EU country.

These are my requirements:

I want to mark each device on the map with a traffic light status indicator - so basically a dot representing each device, it'll be red/orange/green depending on its state.
The states should be updated relatively close to real time, maybe once every 10 seconds if that's reasonable? I'll be pulling state information from a queue that pushes real-time information.
Only details of devices whose states change will be included in each update (Although if required I can push states for all objects visible on the map).
At any one time, ~50k points could need to be displayed, though typically it'll be more like 10k to 15k points. Currently, there are around 30 active jobs with ~175k devices in total (I'd love to show them all, but if it's too much it's not an issue).
There are over 2.5 million points in the database (stored as GPS coordinates) and multiple jobs may be running at the same time, The maps should be able to switch between jobs, but I can limit the number of clients with access to the map, It'll mostly be visible on a display in our team area.
Detailed zoom is not currently a requirement, it'll be a full-screen map on a large display.

So far I'm looking at GeoServer as the server with OpenStreetMap, I think I need OpenLayers too. I've installed postgres and postgis. I will be coding my own service in Java, probably with Spring, but I haven't decided on the front end yet. I'm considering AngularJs and Vaadin at the moment.

So, how realistic is my proposal and what else do I need to consider?


texts:
One of our ArcSDE versions (esri v10.1 and oracle 11g) in SDE.Versions table  shows a 'status' field value of 33554433 instead of the usual 0, 1, 2 for private, public, protected. I understand this is one the SDE workflow status that shows the actual current status of the version in the version's editing processes as discussed here : https://community.esri.com/thread/47331
We have reconciled and posted this version successfully but the status does not seem to change!  Also running the query to show the compression blocking version, shows this as version as the one blocking the compression.
I understand editing the SDE versions table is not a good practice, but I don't see much documentation on how to change the status.  Any advise would be welcome.


texts:
I am not very experienced with the ModelBuilder and am new to ArcGIS Desktop.

I want to use the ModelBuilder to select data from a polygon feature dataset in a first step. The dataset represents the population density of all the districts in a country and I want to only select those, that have a population density of 30% or above.

How can I achieve this with the ModelBuilder?

I am using ArcGIS Desktop 10.4.


texts:
I am new in Google Earth Engine and need help to make the snow line on a merged product of MODIS for aqua and terra. So far, I made the function, union, threshold for NDSI to be greater than 0.4 and a mask, but I need to define the snow line and I don't know how to.
I enclose my code: https://code.earthengine.google.com/fca033adb8e5f3e39e067d1c14d62c69


texts: Splitting polygon from one shapefiles with another shapefiles as cutter using ArcMap

texts:
I want to clip the polygon area. Image below:

using this polygon (other shapefile) image below:

then It will be like this:

I have tried using a clip from the toolbox. but nothing changed.
My plan after that is to add the attribute from the second shp to the first shp.


texts: How to calculate buffer surface that falls within another layer's borders

texts:
The image below shows what two selected layers. One contains a series of borders, and the other a series of buffers.

I am trying to obtain the buffered surface within each of the borders in the other layer, preferably using ftools.


texts: Pipes do not connect to each other - geometric network

texts:
I have created a geometric network in ArcGIS for Desktop. Some of the pipes are not logically connecting to each other even though they appear to geometrically intersect.

Can anyone tell me how to solve this?

For example, I have used "Find connected" at the flag in the figure below. You can see that the pipe is geometrically connected to the intersecting pipe.

I have zoomed into the intersection as much as possible and they seem to be connected.

The figure below show the configuration of the geometric network and the feature dataset.


texts:
I am facing an issue I found in QGIS 2.12.

My layers are WFS layers.

I use QGIS in a Windows tablet to edit some data in the street. After working for some hours in the street and editing a bunch of elements, I toogle the Edit Tool (Switch off) so data is saved in source(client) to postgis database in server (via WFS).

It usually has worked perfectly but this time it did not save the data in the source.  I don't know if the tablet in that momment had a internet coverage problem but the data was not saved in source.

In QGIS, on the other hand, I can see the changes as saved (Edition is switched off) and if I edit another element and save the changes, it is not saved either, and now I have internet connection for sure.

If I close the project and open it again, I am pretty sure I will load everything ok (without the changes of course) and everything will work again, but I don't want to lose the work of a whole day in the street. As said before, all that data is in QGIS project (I haven't close the project).

Is there a way to force all the data in the layer to try to save it to source, even if QGIS now is in a state where it thinks that there is no new changes??


texts: Create x,y points with a radius that best fill a polygon

texts: I am looking to best fill a polygon with non-overlapping circles of a specified radius and then extract the x,y points from a table.  Is there a way to do this in MapInfo without specifically picking the x,y coordinates?


texts: Difference in Bit-depth in Satellite Images

texts: I have DMSP/OLS and VIIRS satellite images. For the 1st image the bit depth is 8, for the 2nd image is 32. I want to make intercalibration between those images. Will it cause a problem due to difference in the bit of the images?


texts: How to get path cost in A* shortest path finder of GeoTools?

texts:
I am computing a shortest path with A* algorithm of GeoTools:

AStarShortestPathFinder pf = new AStarShortestPathFinder(
                                     myGraph, nodeFrom, nodeTo, myAstarfun);
pf.calculate();
Path p = pf.getPath();

How could I then retrieve the cost of the computed path p ?

With DijkstraShortestPathFinder, this can be obtained easily with a pf.getCost(nodeTo) but no similar method seems available with A*. I did not find information in the documentation to do that.


texts:
I am trying to add PostGIS functionality to a database that I created in PostgreSQL using the following command: psql -d <db_name> -f /location/to/this/file/postgis.sql

However, I receive a bunch of errors relating to 'geometry does not exist', but on line 1 of these error messages, it says ERROR: could not access file “$libdir/postgis-1.5”.

I Google'd and checked StackOverflow, and it seems like this error is related to dependencies or scripts not being installed. I installed PostGIS using yum install postgis84. I did not see any dependencies except for proj. Based on these instructions, it seems like I will need Proj4, GEOS, LibXML2, GDAL, JSON-C, etc.

Is there a repository that I can connect to that will let me install PostGIS with all the crucial dependencies rather than installing them one by one? I am running a Red Hat Enterprise Linux 6.5 server on Amazon EC2.


texts:
Novice with ArcGis, I would like to convert my Polyline dataset to a Raster with only 0s and 1s

To this end, I used the Gis command "Polyline to Raster" and set the cell size small enough, so that the lines features might be captured.

However, the Raster Attribute table that I got contained only positively coded cells (i.e. those crossed by a line and whose values correspond to the ID of the line.) while I would have liked to get an attribute table with all the cells.

Thus, I thought about using the "Raster calculator" tool and the following command :

Con(IsNull(raster), 0, 1)

Unfortunately, this didn't work. I got a Raster with an attribute table containing 0 and 1 admittedly, but with only two values now.

Thus, I would like to ask for your help in order to be able to get a Raster with an Attribute table containing all the cells, coded as follows: 0 when no line crosses the cell, and 1 otherwise.


texts: Can't get Tilestache working on Apache (mod_wsgi)

texts:
Trying to get a tileserver running on Ubuntu 14.04 with TileStache and Apache.

No problem running the TileStache development server (tilestache-server.py) or the uWSGI example from the API docs:

uwsgi --http :8080 --eval 'import TileStache; application = TileStache.WSGITileServer("/etc/TileStache/tilestache.cfg")'

However, I can't get it to serve tiles with mod_wsgi through Apache (2.4.7). Followed the instructions here: https://github.com/TileStache/TileStache/issues/67:

Created wsgi file (/var/www/tilestache/tilestache.wsgi, owned by Apache with 0750 permissions) containing:

import os, TileStache
application = TileStache.WSGITileServer('/etc/TileStache/tilestache.cfg')

Following those instructions, added:

WSGIScriptAlias /tilestache /var/www/tilestache/tilestache.wsgi

To my Apache config. Navigating to that alias on any of the defined virtualhosts returns a 403 - Forbidden.

I've also tried replacing that line with the following, but returns the same 403:

WSGIDaemonProcess tilestache
WSGIProcessGroup tilestache
WSGISocketPrefix /var/run/wsgi
WSGIScriptAlias /tilestache /var/www/tilestache/tilestache.wsgi

I've also tried moving those lines inside of a virtualhost, with the same result.

Note: I've changed the tilestache.cfg to use the test/simple cache (not actually caching anything).


texts: Export Rasters from geodatabase containing specific words in title

texts:
I have a file geodatabase containing +200 rasters. I want to export a selection of these based on words in the middle of the file name. Example file names below:

Bat_day_male

Bat_night_male

Frog_day_male

Frog_night_male

Cat_day_male

Cat_night_male

How do I select only the files containing "day" to export to another filegeodatabase? I would like to do this in model builder using an iterator if possible.


texts: Issue running GDAL_Translate on Redhat PPC64

texts:
I try to run the gdal_translate command to extract and convert the MODIS hdf file into its respective band raster tif on a Redhat PPC64 environment and I encountered following error:

with my command:

gdal_translate -sds MOD09A1.A2018273.h27v07.006.2018290185611.hdf MOD09A1.A2018273.h27v07.006.2018290185611.tif

I try to run this command in other environments, it is  working fine. Therefore, I would like to know is there a workaround using python methods that I can look into to extract band raster from MODIS hdf file?


texts: How to retrieve a single polygon from a multipolygon layer in GeoServer

texts:
I have a multipolygon layer in GeoServer. The layer represents states within a country. I want to enable a user to search for a specific state, and then only the selected state to be displayed on the base map.

The problem is two-fold: First, how to retrieve the names of polygons in the layer to enable a user to select, and second, how to display only the selected polygon on the base map.

I'm using OpenLayers for presentation. In addition, I'm using GeoServer and PostGIS. Any pointers on how to proceed?

EDIT: From some research, I'm guessing I'll have to use GeoServer filters (CQL and ECQL) as outlined here. However, I'm still a newbie at this, and any advice will be appreciated.


texts:
I try to draw a graph in map with R

code :

df<-data.frame("from" = c("Lyon", "Toulouse", "Paris", "Marseille"),
               "to"= c("Paris", "Paris", "Marseille", "Toulouse"))
meta <- data.frame("name"=c("Lyon", "Toulouse", "Paris", "Marseille"),
                   "lon"=c(-4.850000, 1.444209, 2.352222, 5.36978),
                   "lat"=c(45.750000, 43.604652, 48.856614, 43.296482))

g <- graph.data.frame(df, directed=FALSE, vertices=meta)
lo <- layout.norm(as.matrix(meta[,2:3]))
plot.igraph(g, layout=lo,vertex.size = 60,
            vertex.color="red",
            vertex.frame.color= "white",
            vertex.label.color = "white",
            vertex.label.family = "sans",
            edge.width=2,
            edge.color="black")

but I do not find a way that allows me to display it on the map using leaflet, someone has an idea please?


texts: Transferring data from PostGIS to Oracle Spatial

texts:
I would like transfer spatial data from PostGIS to Oracle Spatial. I want to use Python's modules psycopg2 and cx_oracle. Transfer Oracle > PostGIS is OK, but in opposite version I have problem with a geometry column. The important part is:

cursorP.execute("select fid, kod, st_asBinary(geom) as geom from ruian.vusc_g")
while 1:
    a = cursorP.fetchmany()
    if not a:
        break
    val = list()
    for row in a:
        dictionary = dict()
        dictionary['1']=row[0]
        dictionary['2']=row[1]
        dictionary['3']=cx_oracle.Binary(row[2])
        val.append(dictionary)
    orcl_inserttext = '''INSERT INTO test_Postgres VALUES (:1, :2, SDO_UTIL.FROM_WKBGEOMETRY(:3))'''
    cursorO.executemany(orcl_inserttext,val)

Wihout geometry it's OK. But in this case I receive Error ORA-01461: can bind a LONG value only for insert into a LONG column. Psycopg2 gives me WKB like buffer, but i don't know, how to use it in a cx_oracle insert.

In my opinion geometry is to large for this syntaxe inserting of records.

Do you have any idea how do it differently?


texts:
I can't find any information about how to use the analyst functionality of Open Trip Planner. The official resources only go on about point-sets and not how to use the framework.

I also looked here: https://opentripplanner.readthedocs.io/en/latest/Analyst/#pointsets
, but it says this documentation is deprecated.

Could someone help me to get the analyst api running?


texts:
I am working on normalizing some rasters so that data taken on different days and calibrations can be seen in the the same relative view.

in gdal_translate the -scale option looks like this:

  -scale [src_min src_max [dst_min dst_max]]:

so if I use the format

  -scale src_min src_max 0 1

I get the result I am looking for now.
However I would like to use the option to not provide the src values and just use the input.

I am wondering if anyone could help me with the proper syntax for this, or if my best option is to just read out the src values before and then apply them in the script.


texts: Deriving slope angle value from DEM using ERDAS Imagine

texts:
I am trying to calculate the steepness of each pixels (30m) from the DEM of 2m resolution. I have made use of ERDAS Imagine - Slope function and returns a contour-like pixel image. The pixel value between contour is null in value. I have also tried to resample the slope image to 30m pixel but i don't know why many pixels return null value even it is overlapped with lots of contour-like pixels.

Is there any method to calculate the steepness of each pixels in the DEM?

What i want to do is to make use of DEM to filter the flat area in the image classification.


texts: Determine minimum number of polygons from shapefile to cover an area of interest

texts:
I have a large number shapefiles representing areas of interest for an analysis that will be conducted using various sources of satellite imagery (IKONOS, RapidEye, etc.). Unfortunately, the imagery doesn't use a pathrow system like Landsat for example, so the extents vary greatly.

I have shapefiles clipped to each AOI representing the extents of different imagery acquisitions, all of which have already been deemed acceptable. Some of these shapefiles have 500 or more polygons.

I need to find an approach, preferably one that can be automated (Python and ArcInfo 10 preferably, FOSS would be acceptable as well) to determine the fewest number of polygons to cover each of my areas of interest.


texts:
I am trying to see if two "circles" intersect. The data I have is  two points and two radii (in meters). I am using srid 4326.  My main issue (i think) is getting the radius in meters.

SELECT (ST_Buffer('POINT (-104.949517 39.740264)',300) &&
          ST_Buffer('POINT (-104.999924 35.000089)',300))

My questions are

Is there a more clever way to do this?
Or how do I convert the radius from dd to meters to get correct buffers?


texts: How to add item to the ArcGIS Portal with the portal API?

texts:
I'm trying to add/update an item back to the ArcGIS Portal using the portal API, in much the same way as you would when doing 'Configure Application'. I am trying to use the update item call to send an item back to the server via POST.

In Chrome I am getting the error "Origin localhost:53657 is not allowed by Access-Control-Allow-Origin". In IE I get the error "Access Denied" and FF nothing happens.

I understand about cross origin, but I thought the ESRI server must support it otherwise how would we make this post? I can't use jsonp since this does not support post

Edit: If I try to do AddItem I get the error:
XMLHttpRequest cannot load https://yyy.maps.arcgis.com/sharing/content/users/yyy_wilton_d/addItem. Request header field X-Requested-With is not allowed by Access-Control-Allow-Headers.

Can anyone see what I am doing wrong? Does anybody have a sample? Or how can I check that the ESRI server is set up to support CORS?

esri.config.defaults.io.corsEnabledServers.push("https://yyy.maps.arcgis.com/");
var apiID  = 07d5de8305cf47c99057ca1336e2yyy
var xhrArgs = {
    url: 'https://yyy.maps.arcgis.com/rest/content/users/yyy_wilton_d/items/' + apiID  + '/update',
    token: this.token,
    overwrite: true,
    content: cont,
    text: dojo.toJson(cont),
    handleAs: "json",
    load: function (data) {
        debugger;
    },
    error: function (error) {
        debugger;
    }
};
dojo.xhrPost(xhrArgs);


texts:
I have problems using the function TIN (linear) interpolation in QGIS. The program crashes when larger areas are interpolated.

Map properties:

vector map,
Geometry typ = Point,
features in this layer: 523851,
EPSG:32633 - WGS 84 / UTM zone 33N,
files = .shp

When I divide the map into more layers (like 80 000 points each) and then interpolate them together (by adding more layers in the interpolation menu), it works, except some areas still makes the program crash.

Any idea how to solve this? Or have someone had/have the same problem?

Maybe the problem is: Can´t handle large files or the algorithm have problems triangulate between some points patterns.


texts: Are there newer routing algorithms (than Dijkstra, A*) in GIS databases?

texts:
There are works like Reach for A* from Microsoft researchers and Highway Hierarchies by Sanders and Schtolz (if I spell the name correctly) from Karlsruhe Uni. Both of them reduce the calculations order a lot, and speed up thousand times on large graphs (see the results in the linked documents). The latter work led to Open Source Routing Machine, which unfortunatly isn't popular enough and not adapted (I couldn't compile it although tried hard).

At the same time, the dbs that I tried, Spatialite and PgRouting, according to their docs, offer just Dijkstra and A* algorithms. I've not even seen bi-directional search mentioned, which saves calculation time twice in my experience.

Are there better algorithms for databases or other applications?


texts:
I want to make a legend on the layer. I heard that in ArcGIS there is a legend wizard, where you can make a legend, but it is in ArcMap, not ArcGlobe.
Is there any tool similar in ArcGlobe?


texts: ArcGIS add attribute table to PointGeometry using arcpy

texts:
I'm trying to add an attribute value to a new PointGeometry before I write it into the new file. How do I do that with arcpy?
Let's say in this example I want to write the ID of the line feature into the new point feature:
points = []
for row in arcpy.da.SearchCursor('lineFeature', ("ID", "SHAPE@")):
    length = int(row[1].length)
    for i in xrange(0, length, 10):
        point = row[1].positionAlongLine(i)
        # should write the ID here
        # something like
        # point.setValue('sourceID', row[0])
        points.append(point)
arcpy.CopyFeatures_management(points, 'pointsFeature')

The command point.setValue('sourceID', row[0]) point. doesn't work, it returns:

Runtime error
Traceback (most recent call last):
File "<string>", line 8, in <module>
AttributeError: 'PointGeometry' object has no attribute 'setValue'


texts:
I have a QGIS layer with points and the column 'name'. This column contains names like "AS 10450 to 1045" or "BAS 1045B100 to AS 1457745".
Now I want to create another column with all values from 'name' but only the part of the string after the 'to', the 'to' should included in the new column.
Some rows does not contain the string 'to', they have values like "AS10" or "SAB145500 14552AS".
I tried the fieldcalculator with different expressions:
"name" LIKE '%to%'

==> I only get the count of the value 'to'
substr("name","name" like '%to%')

==> All values were filled in the new column
regexp_substr("name",'to')

==> I only get the value 'to' in the new column
Where is my fault?


texts: Google Earth won't display a path from a QGIS KML file

texts:
A friend created a KML file(s) for a special bike ride project.  He included a path.  Several operators were using this file in Google Earth.  One however could see everything in the file, mile markers, hospitals, rest stops, etc. but NOT the path.

It appears that this behavior is limited to one specific laptop.  Four others did not have this problem.


texts:
In QGIS 3.0, I am using the new clustering style.

I want to display labels based on cluster size.

What I tried in expression field:

 case
 when @cluster_size>5 then  "label_name"
 end


texts: How to publish metadata in RDF/DCAT endpoint in Geonetwork using custom profile

texts:
I wonder if I can use Geonetwork as a source for harvesting metadata in DCAT-AP-SE (a Swedish profile based on the European format)? I guess there is nothing out of the box, but is it possible to write your own profile plugin, or has anyone already done that and is willing to share? I will use it to publish open data from a Swedish municipaltity.
I've been searching, but the best leads I have found is the Geonetwork documentation on the RDF endpoint that says "in progress" and only gives a brief example of usage, nothing about profiles. I have also found that there seems to have existed a plugin for DCAT-AP, but not for version 2.


texts:
I have a list of vector layers loaded in the canvas of qgis.

  `layers = self.iface.mapCanvas().layers()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.featurecombo.addItem( layer.name(), layer )`

The above snippets of code can be used to display all the layers in the popdown box of the plugin. But it wont load all layers instead it gives results for one layer.(for example, if I have 3 layers loaded and I want a feature count for layers, it will display all 3 layers in the popdown down box but the above snippet will show the feature count  for one layer It would be more helpful if i can get a snippet to load all layers and display the result for all layers at a stretch.


texts: How to install an in house package from local directory in QGIS?

texts:
I need to install a private package which is in my local directory to the QGIS Python environment. I can install it to my conda environmnet by defining the local directory path as channel. How do I install in QGIS Python environment?

Note - I am creating a QGIS plugin which uses services from the in-house package.

The answer to Installing 3rd party python libraries for QGIS on Windows does not work for me as I am not installing Python packages.


texts: I'm having a hard time finding a DEM for the Bahamas. North Bimini specifically. Does anyone have advice? Or a DEM?