texts: Creating new feature class in ArcPy gives error 000210?

texts:
I am trying to create a new attribute table and I keep getting the same error.

import os
import arcpy
import math
folderpath = 'C:\Users\Michaelf\Desktop\GEOG M173'

arcpy.env.workspace = folderpath
arcpy.env.overwriteOutput = True

input_shp = folderpath + '\lower48_county_2012_election.shp'
out_shp = folderpath + '\new_48_county.shp'

arcpy.CopyFeatures_management (input_shp, out_shp)

my error is:

ExecuteError: ERROR 000210: Cannot create output C:\Users\Michaelf\Desktop\GEOG M173new_48_county.shp
Failed to execute (CopyFeatures).


texts: Batch export of layers as individual PNG files in QGIS 3.4?

texts:
I have hundreds of layers. I'd like to be able export each one as a PNG. I can do this by selecting one, and going to project, export map to image, calculate from layer -> select layer -> save. Is there any way I can export all of them as seperate PNG files as a batch process?

I have made progress with the python code below. Unfortunately they are all just white 5k png files.

# create image
img = QImage(QSize(800, 800), QImage.Format_ARGB32_Premultiplied)

# set background color
color = QColor(255, 255, 255, 255)
img.fill(color.rgba())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

# create map settings
ms = QgsMapSettings()
ms.setBackgroundColor(color)

# set layers to render
#layer = QgsProject.instance().mapLayersByName('stream_order')
#ms.setLayers([layer[0]])

#get all map layers:

layers = QgsProject.instance().mapLayers().values()

#loop through layers

for layer in layers:

    # set extent
    rect = QgsRectangle(ms.fullExtent())
    rect.scale(1.1)
    ms.setExtent(rect)

    # set ouptut size
    ms.setOutputSize(img.size())

    # setup qgis map renderer
    render = QgsMapRendererCustomPainterJob(ms, p)
    render.start()
    render.waitForFinished()
    p.end()

    # save the image
    img.save("C:/imageexport/"+layer.name()+".png")

    print('image saved')

replacing rect = QgsRectangle(ms.fullExtent()) with rect = QgsRectangle(layer.extent()) seemed like the logical solution if the issue was the extent being too big, but this just causes QGIS to crash instantly.

I tried the following code and just got a completely black image.

#image_location = os.path.join(QgsProject.instance().homePath(), "render.png")

vlayer = iface.activeLayer()
#vlayer = QgsProject.instance().mapLayersByName("id_")[0]
options = QgsMapSettings()
options.setLayers([vlayer])
options.setBackgroundColor(QColor(255, 255, 255))
options.setOutputSize(QSize(2000, 2000))
options.setExtent(vlayer.extent())

render = QgsMapRendererParallelJob(options)

def finished():
    img = render.renderedImage()
    img.save("C:/Users/smallfing.FINGBIAL/Documents/SHEFFIELD/Road PNG Export/render.png","png")
    #img.save(image_location, "png")
    print("saved")

render.finished.connect(finished)

render.start()


texts:
I'm trying to delete a virtual field from a layer attribute table. The layer source is a Spatialite database so I know that it's not straightforward to remove permanent fields, but as I understand it virtual fields are stored in the project file, not the database.
I followed the answer to this question, but when I try to edit the attribute table the 'Delete field' button remains greyed out:

UPDATE: I've now discovered that the above happens if the layer is added using Layer->Add Layer->Add Vector Layer, but not if using ...->Add Spatialite Layer, in which case a virtual field CAN be deleted using the appropriate button. I'm a relative beginner, and I hadn't realised it was possible to add Spatialite layers in 2 different ways.
However I'm still not clear if this is a bug, or actually correct behaviour, as I don't understand the difference between the two methods of adding layers (though this is probably a separate question).


texts:
Using QGIS 3.4, I am trying to merge multiple DEMs downloaded from the National Map Viewer that are in ArcGRID format (Arc/Info Binary Grid). The output will be a large file so I can't use any of the formats limited to 2 or 4GB in size. I would like to save it as a BigTIFF but I don't know how to distinguish that from a TIFF in the merge raster dialogue box. I tried to save it as TIFF but it failed as file was over 4GB.

How do I save merged rasters as BigTIFF or, if this is not possible, how do I merge rasters to produce the smallest possible file size?


texts:
My work does true ortho rectification.  We are looking to transition into a more open source company.  One of the processes we are looking to "open source" is true ortho.  Does anyone know of a software package (open source) that can perform true orthorectification.  Ideally, it would be a complete package, but it could be a series of libraries, tools, plugins, etc.

Thanks...

Reclarification:

A true ortho is similar, but different to a regular ortho.  With regular orthos you will have building lean.  A true ortho is a process to remove the lean, and "stand-up" the buildings.  The following is an explanation: http://www.esri.com/news/arcuser/1001/standup.html

Because of this, I do not believe this is a duplicated question.


texts: ArcSDE & Oracle & Database Model -- Best diagram the system

texts:
Having Oracle dbms with ArcSDE and ArgGIS Pipeline Database Model (APDM), I would like to best describe the system with a diagram.  In MS SQL each db is autonomous, I think, meaning ArcSDE(using dbo schema) has the SDE business tables in each db and further using a database model, each database has the database model implemented within??

For Oracle (oracle has a schemas instead of db's??), the Oracle instance that has the sde business tables and the database model..assigns users to participate in those tables as a sort of view (schema??) and sde manages those schemas??

My question:  I would like to draw a picture showing how sde and the database model relate as well as how they sit in relation to the rdms on Oracle and on MS SQL Server.  Thanks.


texts: Does Google Maps follow W3C guidelines for accessibility

texts: Does google maps conform with the w3c accessibility recommendations(A,AA,AAA)?


texts: Using GDAL command line to copy projections

texts:
I have two images, a georeferenced Tiff and a second Tiff without a reference. The second image is a copy of the first, but it lost its reference. Now, I want to copy the reference from the first Tiff to the second Tiff using GDAL using the command line.
I had a look at several gdal utilities like gdalwarp, gdal_translate etc. but it seems that they only work on one image at a time and cannot "copy" a projection.
How can I do this?
I am using ImageJ for image analysis (segmentation) on 600 Tiff-images, all formerly geotiffs. However, during import ImageJ sets the xmin,ymin to (0,0) and looses all spatial information. Now, I try to use GDAL within an ImageJ macro to reload the spatial information so that the segmentation output also has a proper spatial reference


texts: Use Round in ArcMap 10.x Field Calculator on field with both numeric and NULL values

texts:
With older versions of ArcMap (9.3, 9.2, etc) I was always able to round numbers down to 2 decimal places using field calculator even if that field had both numeric values and a bunch of NULL values.  It used to just ignore the NULL values.  With ArcGIS 10.x this is no longer the case.

round([UPSTREAMINVERT],2)

Rounding a field that has any NULL values leads to error 999999, error executing function - invalid use of Null: 'round'

I realize it works just fine when you remove NULL values from the selection.  However, I'd prefer not to have to do this for every field I'd like to have rounded.

Has anyone else experienced this and found a work-around?  I've tried to do an if-then statement in field calculator to only return the rounded value if the value is not NULL.  This gives the same error message.

Thanks.


texts: PostGIS installation on Ubuntu 10.12 (Amazon Services)

texts:
On Ubuntu 10.12 (AMI on Amazon), I use aptitude to install Postgresql and PostGis.
On terminal running psql interactive prompt, I type CREATE EXTENSION postgis and

ERROR:  could not open extension control

  file "/usr/share/postgresql/9.1/extension/postgis.control": No such
  file or directory

Is there something I forgot to install?


texts: Application for online overlay images on Google Maps

texts:
I want to develop javascript app which will take any image user uploads, user will click on 3-4 points on the image (3-4 cities, for example NY, Chicago, LA and Miami) and image will be overlayed on Google Maps. I need some expert to tell me if that's possible (something like online georeferencing)?

Thanks.


texts: What are the skills a GIS Expert needs to have?

texts:
The last years I have wondered many times what are the skills a GIS Expert needs to have. I have studied Geology and upon that I did my master in Geoinformatics. Although right now is clear for me what a GIS Expert does, I still find my self wondering how far a GIS guy needs to "go" regarding learning new libraries or frameworks. Perhaps my question can be rephrased to:
Is a GIS Expert a Software Developer? Or should he be?

Its obvious that GIS Expert has to know a lot of different things such

One or two GIS packages
Spatial DBs
Scripting (Python, PHP)
Web Services (WMS etc.)
Web Development Languages/Libraries (for WEBGis)

Actually the above are my skills. It takes lots of time to master them and feel secure using them but I realize that more and more new technologies come in the game.

For example I know how to use OpenLayers API or googlemaps API together with EXTJs and GeoExt. I did several projects on them and I spent considerably amount of time to become good at them.
Now I see other projects using technologies such as GeoDjango or JQuery UI, which I have never used before.

Similarly with PHP. I've been using it for writing scripts to create workflows. Should I go the next level and learn a framework such as Symfony?

Where does the field of a GIS Expert stops and the field of Software Developer starts?

I hope I illustrated my questions clear enough.


texts: To filter the wrong data of point positioning

texts:
I have collected GPGGA data by Novatel OEMV receiver in static single point mode. I have so much data (GPGGA format) and I want to get the static point axis by averaging all data. But some of the data have big deviations so I need to filter the wrong data or data with big deviations before averaging all the data.

In all, how could I filter invalid data? Is there any C arithmetic program or Matlab arithmetic program?


texts: Is there anything like OpenLayers.Filter implemented yet in Openlayers 3? I'm trying to implement a DWITHIN filter using a ol.geom.Polygon and I am searching for API support, which I couldn't find.


texts: GNSS raw measurements on Android: Does real time/post processing Differential techniques possible?

texts:
In the recent months we all saw kind of breakthrough in Android Smartphone GNSS receivers: GNSS raw measurements are now accessible, and this is significant since those are the most common form of widely used GNSS receivers.

https://developer.android.com/guide/topics/sensors/gnss

https://www.gsa.europa.eu/system/files/reports/gnss_raw_measurement_web_0.pdf

After reading the materials, there is one important thing that is still unclear for me: We all know that really precise location can be obtained only by using various differential techniques.

My question is about whether or not the raw data capturing relevant fields used for DGPS (I'm not so familiar with what fields exactly in the GNSS raw data usually used for DGPS):

Does those Android devices can receive real time DGPS messages from existing sources?

And if not, is it possible to construct such service?

For example a simple solution: Place a transportable DGPS reference station near the users, and transmit correction messages as a web service/radio to dedicated client application on Android?

Also, does the raw measurement data available on Android, can be stored for later post processing (offline) DGPS correction?

This is relevant in case different raw measurements used for real time/offline, and maybe the available fields in the raw data on Android are only those used in offline.


texts: Feature "Add to wishlist" in Info Window in Google Maps API 3?

texts: I have Google Map and a lot of markers. Is there possibility to add feature when user clicks on marker and Info Window shows up, that user has option to add that location to wishlist and associate that to either a user ID/cookie?


texts: Define map palette using variable from FeatureCollection in Google Earth Engine

texts:
I'm trying to plot the ESA CCI land cover map using the palette that they define in their .csv download. I managed to import this palette into GEE by saving it as a shapefile with a dummy geometry (using R), which I'm hoping I've made publicly available here.

I'm testing it out by changing the palette in GEE's Image Visualization tutorial (for a different land cover dataset, but I'll stick with it for simplicity). So, first, I load the MODIS land cover:

// Load 2012 MODIS land cover and select the IGBP classification.
var cover = ee.Image('MODIS/051/MCD12Q1/2012_01_01')
  .select('Land_Cover_Type_1');

Then I define my alternative palette:

// Load alternative colour palette
var cciPalette =
  ee.List(
  ee.FeatureCollection('users/rasenior/cci_landcover_legend')
 .aggregate_array('hex'));

// Remove '#' symbol and make lower case
var replace = function(name){
  return (ee.String(name).replace('#', '').toLowerCase());
};
var cciPalette2 = cciPalette.map(replace);
print(cciPalette2);

The console tells me that cciPalette2 is a list, and everything looks good. It looks just like the igbpPalette from the tutorial (a longer list because CCI have more land cover categories). But when I sub cciPalette2 in for igbpPalette:

Map.setCenter(-99.229, 40.413, 5);
Map.addLayer(cover,
             {min: 0, max: 17, palette: cciPalette2},
             'CCI classification');

I get the error:

CCI classification: Layer error: Image.visualize: Expected a string or list of strings for field 'palette'.

But cciPalette2 is a list of strings, no? What am I doing wrong? I'm a total GEE and javascript novice so I'm sure I'm missing something simple.

Full code here:
https://code.earthengine.google.com/be1102eda2fb41584e9be3697f42b48e


texts: OpenLayers v5 XYZ source is making 'wrong' requests

texts:
I'm trying to add a tilelayer to a map, but the zoom level and x/y coords are not matching up in the network request for some reason.

Here's my code for the layer:

new TileLayer({
  source: new XYZ({
    url: 'myserver.com/tms/{z}/{x}/{-y}.png'
})

In the browser's network tab, the request looks something like myserver.com/tms/10/585/716.png and the response is an empty png file.

However when I take that url and manually decrease the zoom by 1, so it becomes myserver.com/tms/9/585/716.png I get the correct tile as a response. This is consistent across all zoom levels.

For reference, the OSM source layer works perfectly fine:

new TileLayer({
  source: new OSM(),
})

I'm guessing this is an issue on my part since the server is technically serving the right tiles but other than that I don't know what to do.


texts: Use open layers plugin display as georeferenced image in QGIS

texts:
Could someone tell me if it is possible to save a
plugin open layers display as georeferenced image in QGIS 2.18?


texts: Managing members in ArcGIS Online - Multiple member content transfer then delete

texts: I would like to delete about 15+ members in our organization's ArcGIS Online account. I know I can transfer the content of each member to another member's account and then delete the former. I would like to know if their is a way to do this for multiple members at one time?


texts: QGIS - Calculate Highest Elevation within 1 degree grid

texts: I have downloaded country raster elevation data from Diva-GIS and added 1 degree grid from "View >> Decorations >> Grid" as seen in the image below. Now, I want to know the highest elevation within each grid, how do I get this done?


texts: Randomly Generate Streeview Point without GIS in country

texts:
I'm interested in pulling street view data at random within a country.

I can only figure out how to generate a point within a bounding box. Often it pulls data from nearby nations. I'm using a GIS server right now, but I'm interested in removing it from the workflow since it's a pain.

Anybody know of any javascript mapping API that can randomly generate a point within a complex polygon, such as a country boundary?


texts: How to execute sql lines in postgresql using python？

texts:
I use postgresql 8.4 to route a river network, the following is my code typed in pgAdmin and the result is fine,

select * from driving_distance
('select gid as id,
 start_id::int4 as source, end_id::int4 as target,
 shape_leng::double precision as cost from network',
 10, 1000000, false, false);

I have installed psycopg2 and it successfully connect python and postgresql,

#set up python and postgresql connection
import psycopg2

try:
    conn = psycopg2.connect("dbname = 'routing_template' user = 'postgres' host = 'localhost' password = '****'")
except:
    print 'I am unable to connect the database'

Now I need to directly execute my sql code on the top in pyscripter, how should I change these codes to python code？

I am working with postgresql 8.4, python 2.7.6 under Windows 8.1 x64.


texts: QGIS Print Composer - Display Symbols for Empty Classes?

texts:
I am working on a map of mortality rates for several racial/ethnic groups and I have created custom classes to encompass the rates instead of the QGIS classes created through its various modes. This is to create a standard set of groupings and to minimize confusion created by different ranges being given similar coloring schemes across groups.

The problem I am having lies in the fact that some groups are not represented in all the classes. For instance, a range from 0 to 1200 with classes of 300 (0-300, 300-600, and so on).  One group never goes above 600 while another group never falls below that level.  And the print composer, at least as far as I've been able to identify, will not include those classes which are empty.

To allow for a visual comparison between the maps, however, I would like the legend symbols to be present even if they are not used in the map.  The reader can parse from the colors on the map that some of the categories are not fulfilled and it allows for easier interpretation for a reader comparing the maps side-to-side.

I'm hoping this is a simple issue where I've just missed something.


texts: Setting up QGIS ECW support on Linux Arch / Manjaro?

texts:
How do I setup ECW support for QGIS in Linux?

In windows it's super easy, due to the OSGeo4W installer, but I can't get it running in linux.

I tried a lot of topics from gis.stackexchange like these:

ECW for QGIS 1.8 in Ubuntu 12.04
Unable to install ECW support on lubuntu 14.04

or other blogs like this one:
3

The main problem is, that the AUR package for libecwj2 isn't online anymore and I'm still new to compiling software in linux. So somehow the compilation, which is described in the first two links, always fails.

@SamTux replied but I still can't get it running.

The archgeotux-repository did work so far, but I somehow I got the following error: link
I wasn't able to open qgis anymore and it looks like a lot of more trouble. So I thought this can't be the answer.

I could install the libecwj2 AUR package (YUHUUU!), but got stuck with the recompilation of gdal. I found it also strange why it's not possible to edit PKGBUILD of the gdal but for the gdal-filegdb package?

When I do the recompile of gdal-filegdb the package "libkml" fails: build-log

Another question, when I upgrade in the future to a newer version of gdal, do I always have to edit the PKGBUILD do ensure the ecw support?


texts: How to display graduated symbol in the legend in print composer?

texts: This is the problem: why the graduated symbol is not reported at all in the legend of the print composer? Is it a bug or am I doing something wrong? Here below an image that shows the problem:


texts: QGIS 3.4.5 print composer: doesn't print DIN A0

texts:
I want to find a way to print a map with a raster layer in greyscale (GeoTIFF) and two vector layers in shp-format.

I created a map layout and went to the page properties (right-click on the canvas) to enter the Page Size A0. Exporting the layout to pdf results in a map DIN A4.

Can somebody tell me about the settings I have to change in detail so that I can plot my map?


texts: ArcGIS 10.5 auto-populating field in attribute table

texts:
I am having trouble auto-populating a field in the attribute table.

I would like to populate the grid_code column with sequential letters and numbers - i.e. AAA001, AAA002 etc. I looked at previous threads on auto-populating, but was unable to get it to work.

Image of attribute table attached:


texts: Using annotation layer inside group layer without getting error adding element?

texts:
How do you get an annotation layer working inside a group layer? When I try, I get the error message "error adding element".

I have a project that requires several annotation layer to be organised within a Group Layer.

Tested in ArcMap 9.3.1 and 10.

Steps to reproduce:

Create a new document
Add a Annotation Featureclass to ArcMap
Start Editing
On the Draw Toolbar, change the Active Annotation Target to the Annotation Featureclass
Draw something (e.g. freehand) - This should work!
Create a New Group Layer
Move the Annotation Featureclass inside it.
The Annotation Target has now been lost!
On the Draw Toolbar, change the Active Annotation Target to the Annotation Featureclass
Draw something (e.g. freehand) - This does not work!
You will get the error message "error adding element"

You will now not be able to write to this Featureclass unless you save the document, close ArcMap and re-open it. This is not an acceptable workaround for my project. Is there a workaround, preferably in ArcObjects?

This might relate to an ArcMap 8.3 bug submitted in 2001.


texts: Rasterize in QGIS using point as centroid (not vertex)

texts:
I'm attempting to make the full switch from ArcGIS to QGIS, but there are a few snags I'm hitting along the way, one of which is in the GDAL "Rasterize" function. My process goes like this: I have a raster which I convert to point file using SAGA "Raster values to points" which creates a point at the centroid of each pixel (with XY coordinates as attributes). I then run a function on this dbf outside of QGIS (in R). Then I want to convert the point file back into a raster. I use the GDAL "Rasterize", but this uses the points as vertices, effectively shifting the raster by half a pixel down and to the right. By comparison, ArcGIS views the point as the centroid, and using the "point to raster" returns a raster that lines up with the original. How can I tell QGIS that the point is the centroid?

I have tried the roundabout way of making a VRT, but as the pixels are shifted exactly half a pixel, the resultant raster is shifted an additional half-pixel down and to the right (so now the pixel vertices line up, but are shifted a complete pixel diagonally down and to the right).


texts: OpenGeo Suite 3.1 doesn't work anymore when setting up Client SDK 3.0.2

texts:
We've installed OpenGeo Suite 3.1 on Ubuntu 12.04 LTS and everything works fine.

When trying to set up the OpenGeo Client SDK 3.0.2 (following http://suite.opengeo.org/opengeo-docs/apps/clientsdk.html), however, the OpenGeo Suite doesn't work anymore.
After inserting the SDK in PATH and adding ANT_HOME and JAVA_HOME in /etc/environment, tomcat6 starts up, but doesn't serve anything in the browser anymore.

The severe error in the log looks as follows:

24 Jul 15:35:24 WARN [config.CustomEditorConfigurer] - Passing PropertyEditor instances into CustomEditorConfigurer is deprecated: use PropertyEditorRegistrars or PropertyEditor class names instead. Offending key [org.geotools.util.Version; offending editor instance: org.geoserver.platform.util.VersionPropertyEditor@3f705316

24.07.2013 15:35:24 org.geoserver.platform.GeoServerExtensions checkContext
SCHWERWIEGEND: Extension lookup occured, but ApplicationContext is unset.

24 Jul 15:35:24 ERROR [geoserver.global] -
- GEOSERVER_DATA_DIR: /usr/share/opengeo-suite-data/geoserver_data

Any help would be greatly appreciated. Are there conflicts between the versions? A new SDK doesn't seem to exist yet.

EDIT: The error from the log doesn't seem to have to do anything with our problem, since it is also logged in another working geoserver environment.


texts: Calculate bearing between two decimal GPS coordinates

texts:
I'm building a virtual tour application that uses custom panoramas in conjunction with Google's StreetView service.  With StreetView, you can have 1 to n links at each particular location to move to, and these links are based on compass bearings.

Each panorama has an associated GPS location to six decimal places acquired via iTouchMaps.  My issue is when I calculate the compass bearings from each link location, it returns a bearing that is a tiny (as in 0.00001 degree difference) from true north.

Given a map like

                       4
                       |
                       |
                       |
                       2
                       |
                       |
                       |
                       1

where each number is a location and has the following values:

1- 43.682213, -70.450696
2- 43.682194, -70.450769
4- 43.682179, -70.450837

What I'd like is a method to determine programatically that the bearing from 1 to 2 and 2 to 4 is 70 degrees, and 2 to 1 is 250 degrees.  Using the code from here as a starting point, I wrote the following function

//starting lat/long along with converting degrees to radius
    var endLat = location.lat();  //for rhumb calc
    var endLong = location.lng();

    //loop over response, calculate new headings for links and add link to array
    for(var i=0; i<data.length; i++){
      //this link's lat/long coordinates
      var startLat = data[i].lat;
      var startLong = data[i].lon;

      //get the delta values between start and end coordinates
      var dLong = endLong - startLong;

      //calculate rhumb bearing to
      var dPhi = Math.log(Math.tan(endLat/2+Math.PI/4)/Math.tan(startLat/2+Math.PI/4));
      if (Math.abs(dLong) > Math.PI ) dLong = dLong > 0 ? -(2*Math.PI-dLong) : (2*Math.PI+dLong);
      var bearing = (toDeg(Math.atan(dLong, dPhi)) + 360) % 360;

Using the above function and substituting

location for position 2
data[0] for position 1
data[1] for position 4

I get the bearings

359.9958174081038 degrees to 4 (should be 70 degrees)

and

0.0038961130024404156 degrees to 1 (should be 250)

Why is this function returning such different values than expected?


texts:
I'm composing a print layout with a map on QGIS, and I'm trying to shape the map based on a layer.

I have a layer that sets the limit of the study area, and I need to put on the layout only the things inside the limit. Outside the limit of the study area I want it to be empty, without the satellite image and other informations.

Does anybody knows how to do that?


texts: Problem in Adding line feature in Openlayers vector layer

texts:
I have added the WMS line feature in vector layer. But I am facing some display problem in it as shown below screenshot.Is there anyother way to do this?

But there is no issue in point features.

Here is my code:

 selectionLayer = new OpenLayers.Layer.Vector("Selection", { styleMap:
            new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"]),
    displayInLayerSwitcher: false,
    isBaseLayer: false
});

osMap.addLayers([selectionLayer]);

var wmsFeatInfo = new OpenLayers.Control.WMSGetFeatureInfo({
url: 'http://xxx.xxx.x.xx:8080/geoserver/cite/wms',
    title: 'identify features on click',
    queryVisible: true
});
wmsFeatInfo.infoFormat = 'application/vnd.ogc.gml';
wmsFeatInfo.events.register("getfeatureinfo", this, getWmsFeatAttributes);
osMap.addControl(wmsFeatInfo);
wmsFeatInfo.activate();

function getWmsFeatAttributes(e) {
    selectFeatLst = {};
    if (e.features.length > 0) {
        for (var ii = 0; ii < e.features.length; ii++) {
           var layer_ID = e.features[ii].fid.split('.');
           var featId = layer_ID[1];
           selectFeatLst[featId] = e.features[ii];
        }
    }
 }

 function test(key){
    selectionLayer.addFeatures([selectFeatLst[key]]); //Here line feature is adding
 }

Please suggest me to fix this issue.


texts: How to calculate midpoint given north/south latitude and east/west longitude?

texts:
First, sorry if this questions sound very simple. I am very new to GIS.

Given four values of north/south latitude, east/west longitude how to calculate the center point that is correct for all cases on earth?

                  ----------(north latitude)---------------
                 |                                         |
                 |                                         |
 (west longitude)|                ? (calculate this point) |(east longitude)
                 |                                         |
                 |                                         |
                  ----------(south latitude) --------------

The dateline is an issue if you just use the following formula (N+S)/2 and (E+W)/2 to calculate the center.


texts: How to add my photo in KML with Android?

texts:
I have tried with overlay and Placemark, I store my photo and KML together in one KMZ but when I open with Google Earth,  I only can see the text in the location, not the photo. Here is my code:

[?xml version='1.0' encoding='UTF-8'?]
[kml xmlns='http://www.opengis.net/kml/2.2']
[Document]
[name] Bungay[/name]
[description]Photos.  Walking distance from Home.
[/description]
[PhotoOverlay id="photo1"]
[name]
Garden1.jpg
[/name]
[Icon]
[href]
Garden1.jpg
[/href]
[/Icon]
[Point]
[coordinates]
1.4517828,52.4515304,0
[/coordinates]
[/Point]
[/PhotoOverlay]
[/Document]

[/kml]

And with Placemark :

[?xml version='1.0' encoding='UTF-8'?]
[kml xmlns='http://www.opengis.net/kml/2.2']
[Document]
[name] Bungay[/name]
[description] Photos.  Walking distance from Home.
[/description]
[Placemark]
[Icon]
[href] Garden1.jpg[/href]
[/Icon]
[name] Garden1[/name]
[ExtendedData]
[Data name='Group']
[value] Garden [/value]
[/Data]
[/ExtendedData]
[address] [/address]
[Point]
[coordinates>1.4517828,52.4515304,0
[/coordinates]
[/Point]
[/Placemark]
[/Document]

[/kml]


texts: QGIS: Converting Projection Coordinate Systems to Geographic Coordinate Systems to replace (X,Y) with Latitude and Longitude

texts:
So as the title suggests, I have a shapefile that shows a set of rail lines and I want to display each line independently so I extracted the nodes with Geometry Tools. Unfortunately these coordinates are in PCS and I need them in GCS. Is there an efficient way to do this? Like a plugin? I have 10xxx points and am unfamiliar with Python. (590462.0801,14157127.48) and I need something like (39.508974,-119.814971).

As you can see from the coordinate window, I'm still stuck with a PCS value.  I've successfully changed my CRS but it hasn't had an effect.


texts: National Land Cover Dataset (USA) Data Extraction

texts:
I have downloaded the NLCD 2011. It has different formats of files (i.e. .rrd, .ige). I tried to open the (.rrd) file using ArcGIS 10.5.1 but it couldn't open it.

Can you suggest me how can I view the land cover classes from here?


texts: OpenLayers 3 OSM layers in ExtJS 5 looks stretched

texts:
I've created OL3 OSM layer in ExtJS 5 (previously did it with OL2).
On my desktop everything looks fine , but on the mobile phone the map looks stretched (This is denver for example).

This is the code I used :

Ext.define('ExtOL.view.Map', {
extend: 'Ext.panel.Panel',
region: 'center',
xtype: 'map',
rtl: true,
html: "

", // The map will be drawn inside
stateful: false,
border: true,
listeners: {
    afterrender: function () {

        var extent = ol.proj.transform([34.7193330, 31.8573201, 34.9892725, 32.2086729], 'EPSG:4326', 'EPSG:3857');
        var center = ol.proj.transform([34.8543028, 32.0329965], 'EPSG:4326', 'EPSG:3857');

        layers = [];

        var osm = new ol.layer.Tile({
            source: new ol.source.OSM({})
        });

        layers.push(osm);

        var view = new ol.View({
            projection: 'EPSG:3857',
            center: center,
            zoom: 15,
            minZoom: 0,
            maxZoom: 16
        });

        map = new ol.Map({
            target: 'gis_map',
            renderer: 'canvas',
            layers: layers,
            view: view
        });

        map.getView().fitExtent(extent, map.getSize());

    },
    // The resize handle is necessary to set the map!
    resize: function () {
        var size = [document.getElementById(this.id + "-body").offsetWidth, document.getElementById(this.id + "-body").offsetHeight];
        //console.log(size);
        map.setSize(size);
    }
},
   initComponent: function () {

        this.callParent(arguments);
    }
});

And this is the App init :

Ext.define('ExtOL.Application', {
extend: 'Ext.app.Application',
name: 'ExtOL',
views: ['Map', 'Panel','LayerSwitcher'],
controllers: [],
stores: [],
launch: function () {
    Ext.create('Ext.container.Viewport', {
        layout: 'border',
        rtl: true,
        style: 'border:1px solid;padding:2px;background:#157fcc',
        bodyBorder: true,
        defaults: {
            bodyPadding: 0
        },
        items: [
            {
                xtype: 'map'
            }
          ]
        });
    }
});

Any idea why the map is stretching ?


texts: Handling events that cross more than one route using Linear Referencing of QGIS

texts:
I am working on implementing linear referencing.

How would you deal with the situation below?

Where an event is stretched across two or more "routes". At the moment I am playing around with the LRS plugin in QGIS, but hoping to implement it using postGIS only when I am done.

Another similar situation would be road maintenance, where the same project is to improve the road surface for a section that cross multiple roads. So the starting point of a an event can be on route 2 at X meter and end at route 1 at x meter. The examples I have come across so far is using event layers where the event start and stop on the same route, I find it very limiting.


texts: Arcpy: MXD automatically removing layer whose data source was overwritten

texts:
I have a large script that does many things, but at the end it creates a buffer whose output is an existing feature class that is overwritten.  That feature class is the data source for a layer that is in a mxd that is current an object in memory.  For some reason, that layer from the mxd is being removed as soon as I run the buffer.  I isolated the problem with the code below by saving out two mxds.  When I open the 1st mxd has the layer in the TOC, but in the 2nd mxd it is gone.

sqlWellName = "NAME = " + "'" + WELL_NAME + "'"
SDE_lyr = "SDElyr"
arcpy.MakeFeatureLayer_management(r"C:\GDB_Connect\Editor\MAIN@DEFAULT.sde\WATER.WELLS", SDE_lyr)
arcpy.SelectLayerByAttribute_management(SDE_lyr, "NEW_SELECTION", sqlWellName)
mxd.saveACopy("C:\ProjectData\Map_1.mxd")
arcpy.Buffer_analysis(in_features=SDE_lyr, out_feature_class=r"C:\ProjectData\DATA\AutoMap.gdb\Buffer1Mile", buffer_distance_or_field="1 Miles", line_side="FULL", line_end_type="ROUND", dissolve_option="NONE", dissolve_field="", method="PLANAR")
newmxd.saveACopy("C:\ProjectData\Map_2.mxd")

Is this a bug?


texts: How to calculate the area of each WRB soil class with the 2017 TAXNWRB_250m raster?

texts:
For my research I want to calculate the total area of each WRB soil class on each continent in km2 using the 250m WRB soilgrids raster. However, the total area I calculated is almost 1.3 times as large as the total land surface area on earth.

I downloaded the TAXNWRB_250m_1l (2017 data) raster from here https://files.isric.org/soilgrids/former/2017-03-10/data/ and loaded it in arcgis pro. I used "extract by mask" to create a raster for each continent, using the continents as masks.

I assumed that because the resulution is 250m, I could calculate the total area in km2 by multiplying the number of cells for each soil class with (0.25*0.25). I used the "Count" column in the attribute table for the number of cells per soil class. However, there is a total of 3,275,963,752 cells for all continents together, which gives a total area of 204 million km2.

This is remarkable since the total land surface area  is about 148 million km2. The original raster I downloaded has even more cells (3,280,092,749), so that would even more overestimate the surface area.

Area = number of cells* (cell_size*cell_size) is common practice I believe, so this makes me think it has to do with the cell size I use.

Could it be that this has something to do with the projection, so that the 250m resolution only applies at the equator and becomes smaller at higher lattitudes? The raster only has a GCS (WGS84) and I tried projecting it to both Mollweide and Homolosine, which give 217.7 and 224.3 as cell sizes, respectively, when I look into the properties. Using these as cell size the total area still exceeds the theoretical land surface area of 148 million km2.

Any advice on how to correctly calculate these areas would be most welcome!


texts: PostGIS created fixed dimension raster bigger than geometry

texts:
I'm having trouble creating a raster as a feature map for my algorithm.

My algorithm needs a numpy array of fixed dimension (e.g. 200x200) indicating if there is land (1) or no land (0) at that position.

I have a PostgreSQL database of Heteogeneous Integrated Dataset https://core.ac.uk/display/158460560 which contains a table named 'geographic_features.europe_coastline_polygone'.

I select geometry for my specified location xmin, ymin, xmax, ymax using a intersection query and this seems to be fine.

But now I want to rasterize the geometry. Since I've found that creating an empty raster and using that as template I'm doing that.

ST_MakeEmptyRaster({width}, {height}, {xmin}, {ymax}, {x_width}, {y_width}, 0, 0, 4326)

All of these parameters are calculated in python. x_width and y_width are calculated as x_width = (xmax - xmin) / width.

Then I use ST_asRaster(geom, template_raster)

The loading is done via WKB using the second answer from here Downloading raster data into python from postgis using psycopg2

My problem is that I'm not getting a array of size (width, height) = (20, 20) as expected, but an array of size (12, 12).

I'm assuming that PostGIS only rasters the bounding box of the geometry instead of the specified width and height.
Am I wrong here? Searching the problem only yielded questions about fixing the cell width, which seems to be working fine already for me.


texts: I need to perform Gram-Schmidt pan-sharpening, I don't have access to ArcGIS or ENVI so have been using QGIS. Is there any plug-in or similar capability that is available to perform Gram-Schmidt pan-sharpening in QGIS? If not, is this possible in any other open-source software?


texts: Get messages from erroneous python script run in C# .NET with ArcObjects

texts:
I'm trying to implement some error handling into my C# code for my .NET application. I feel like I'm missing something really obvious because I know there has to be a way to do this, but I just can't figure it out. I have a python script I call out to using gp.Execute(). If there is an error in the python script and it does not complete, I cannot get the message from it, because any call like gp.GetMessages() has to come after the gp.Execute() statement, so it never gets hit if the script encounters an error and execution stops in my C# code. How do you get error messages from your scripts? I have also tried to put gp.GetMessages() in the catch statement but the gp object is considered out of scope there.

My current code:

try {
     //stuff here

     //Run PopulatePatches script
     IGeoProcessor2 gp = new GeoProcessorClass();
     gp.AddToolbox(@"C:\Users\mxxx\Desktop\GIS_Testing\HbtatTesting\HabitatMon.tbx");
     IVariantArray parameters = new VarArrayClass();
     parameters.Add(shapefile);
     gp.Execute("PopulatePatches", parameters, null);
     gp.GetMessages(ref sev) //never hits this line if python script fails
}
catch {
     ViewBag.Message = "Failed!";
     return View();
}

Lastly, I have gotten rid of the try/catch block and tried to just get the script to play through by doing this:

IGeoProcessorResult pResult = gp.Execute("PopulatePatches", parameters, null)
int severity = 2;
string messages = pResult.GetMessages(severity);

//I have also attempted the above line with just straight gp.GetMessages(severity)
//because I think result objects are null if the tool doesn't execute properly

ViewBag.Message = messages;
return View();

It will always stop on the execute call, throwing the general HResult E_FAIL COMException in .NET even though I know the error is a ZeroDivisionError in the script because I purposely put it in there (the code runs perfect in python and my .NET app if I fix the zero division error).


texts: ArcGIS Desktop enforcing Vertical Coordinate System properties for rasters?

texts:
In upgrading from ArcGIS 10.4.1 to 10.5.1 I've noticed that ArcMap seems to handle raster coordinate system info a little differently now.

We have an older set of 10m NED elevation rasters on our data server. In 10.4.1 these files could be re projected from their native GCS NAD 83 to any other common projected system. No errors.

In 10.5.1 the same input raster will not project because Arc interprets the vertical coordinate system (VCS) as "unknown" and the project raster tool hangs.

If I check the reference system in the image properties in 10.4.1 I can't see any lines for the VCS. In 10.5.1 there are VCS properties for the same image and they are set to "unknown" (unknown VCS datum...). If I redefine the input raster and set the VCS on the Z coordinate system tab all works ok in 10.5.1 even if the output of project raster is set to have no VCS!

The issue is, I don't understand why ArcMap suddenly cares if the VCS is set?

The input file is the same in both 10.4.1 and 10.5.1. In 10.4.1 Arc seems to interpret an unknown VCS as "no VCS", in 10.5.1 it seems like I have to explicitly tell it it's no VCS or set a real VCS.

What has changed about VCS handling in ArcMap at 10.5.1?


texts: GDAL couldn't create raster with equal shape

texts:
I am trying to carry out pixel by pixel correlation for each image in a time series vs another 1d array of data. To do that, the arrays coming out of the raster files must be identical in shape. However, when I use the code below, I get different shapes.

def change_raster_extent(extent, in_file, arrays):
    out_tif = in_file.replace('.tif', '_ext.tif')
    ds = gdal.Open(in_file)
    gdal.Translate(out_tif, ds, outputBounds=extent)

    band = ds.GetRasterBand(1)
    # dis_arr = np.reshape(discharge, (21, 0)).T
    x_array = np.array(band.ReadAsArray())
    print x_array.shape
    arrays.append(x_array)
    return out_tif

extent = [-98.05392182  17.40680802 -80.43304074  31.46484375]
arrays = []
for in_file in in_files:
    change_raster_extent(extent, in_file, arrays)

It outputs:

(971, 1217)
(982, 1231)
(985, 1234)
(985, 1235)
(991, 1242)
(995, 1247)
(993, 1245)
(967, 1213)
(992, 1244)
(978, 1226)
(986, 1236)
(972, 1219)

Edit: I realized the possible cause of this issue is because these images are MODIS aqua satellite data and their pixel positions don't match. So even though their geographic extent is same, the pixels don match. I suspect, re-sampling the images could fix the issue. I will try that and see if it fixes it.

How could I make the shape to be identical for the output array?


texts: Do not copy domain with the feature dataset

texts:
Is there a way to copy a feature dataset to another database with all the features but not carry the domain with it?  Version 10.4.1.  The problem is that I only maintain the domains in the Parent database.  If a domain is edited and then different then the child database, then Arc creates new domain in the child with an _1.  I do not want to have to maintain the domains in all the child databases.

I hope there is a work around this problem


texts: Seeking datasets to construct minimal reproducible examples for ArcMap

texts:
I would like to construct a minimal reproducible example to illustrate a question here on GIS SE. Ideally, that example would reference data that's easily accessible to any ArcGIS user, without requiring that they download data from some external site.

Are there example feature classes (and/or data of other types) that I can rely on every user having access to as part of their ArcGIS Desktop (ArcMap) installation?


texts: How to draw a polygon given a distance and bearing in QGIS?

texts: I want to draw a polygon given a length and a bearing, basically a square with length 1m. I can't seem to manage it with the existing tools in QGIS 2.4 (I possibly need a plugin).


texts: Map services disappear in Silverlight and Internet Explorer

texts: I upgraded from Silverlight Viewer 3.0 to 3.2 and upgraded one of my projects.  Now in that project everything disappears at the last 2 zoom levels.  If I switch to dynamic instead of cached it displays.  But I tried recaching and it didn't help.  It will also display fine in Mozilla.  Any idea why IE (10 or 11) are behaving in this manner?


texts:
My maps georeference just fine in QGIS and I save them to jpg with jpgw.  The jpg image loads but the jpgw file does not.  I keep getting the error message "this is an unsupported file", or words to that effect.  I have tried bmp and bmpw, tif and tifw and tiff and tiffw with the exact same result.

I got the same result with Global Mapper but surprisingly, OK Maps loaded both the image file and the world file perfectly.  My GPS receiver tracks accurately with OK maps also.

Please help any way you can.  I really enjoy using QGIS.

Thanks, Jim


texts:
I am experiencing some shifting polygon that I obtain via pixel to lat lon conversion. I have a program to do image processing in opencv that uses pixel coordinates and once I had obtain my contours, I would convert these into polygon by transforming the pixel points in the polygon to lat Lon and visualized on QGIS.

I am stuck now because I keep experiencing a shift in the polygon when plotted on lat Lon epsg 4326 when I plot the polygon on my original georeference image geotiff.

I suspect it may be a precision issue but seeing that the shift is 10-9 pixel in x and y direction which translate to about 0.1 - 0.2 m error on the geotiff file, I don’t think that’s the case.

Process flow:

Geotiff image
gdaltranslate to generate jgw, prj, jpg.
process jpg
Read prj and jgw for georeferencing and transformation between pixel to lat lon and vice versa.
obtain contours
transform polygon
convert from pixel to lat Lon
overlay polygon on geotiff to visualize


texts: Creating a mosaic from multiple Landsat scenes on Qgis

texts:
I wish to produce a large image. As a base use cloud free regions of multiple Landsat images. I wish to keep the bands and be able to work with the generated mosaic just like a "normal" Landsat image.
What ways would you propose to archive a satisfying result?
I am on Linux/Qgis


texts: Bending arrowlines between symbols (points) using MapServer

texts:
Good day everyone. I need help of the community, and asking for it:

There is such a task: to draw a curved line with an arrow at the end between two symbols (points).

Questions:

What is the best way to bend a line defined by two points? (current solution: I find the coordinates of the third point, randomly adding from 0.25 to 0.75 of the difference between the ordinates and abscissas of the points, then I break the line into smaller segments and interpolate it. But since the difference between the ordinates or abscissas can be insignificant, the line can almost do not bend).
How to draw arrows at the ends of lines using MapServer? (current solution: using GEOMTRANSFORM, add a filled-in character 'V' to the end of the line, find the angle between the last two points to rotate the (POLAROFFSET) character so that the arrow points to the center of the character at the end of the line. Main problem: when scaling characters, the arrow shifts not along the line, but along the angle that is given, and since the line is a curve - the angle changes and the arrow moves off the line).
Add. clustering question: Is it possible, when clustering the points of the beginning and end of lines between themselves, to replace the lines with the line between the clusters (having new coordinates (between points in the cluster))?

MapServer 7.7-dev / 7.6.

Approximate image:


texts: How to set project CRS with plugin (sqlite layer)?

texts:
I'm loading a sqlite db into the current project of qgis with a plugin. The standard CRS is EPSG:4326, however i want the canvas to be set at EPSG:3857.
When running one of the 2 codes in the console, it works correctly.

After running this code:

# Settign the CRS of the data set option 1
print('setting CRS to EPSG 3857')
QgsProject.instance().setCrs(QgsCoordinateReferenceSystem(3857))

I can see a very quick change of the CRS code in the right lower corner of qgis. So i gues it changes to EPSG:3857 but it changes back to EPSG:4326. But when i reload the layer again (so by now there are 2 the same sqlite layers in the toc) it changes to the correct EPS:3857 CRS.

The code below from:

# Settign the CRS of the data set option 1
canvas = self.iface.mapCanvas()
selectedcrs = "EPSG:3857"
target_crs = QgsCoordinateReferenceSystem()
target_crs.createFromUserInput(selectedcrs)
canvas.setDestinationCrs(target_crs)

With this code the CRS is correct, EPSG:3857, but in the lower right corner the label states the wrong CRS code. EPSG:4326 in stead of EPSG:3857. The data is for sure EPSG:3857.

Main question:
how do i load the sqlte db into the current project with the correct CRS?


texts: Creating spherical texture map from QGIS World map?

texts:
Right now I am building an application in OpenGL. The application maps a spherical Earth texture map to a polygon approximating a sphere. When I use any of the sphere textures such as from http://planetpixelemporium.com, the map fits neatly and latitude and longitude values returned are accurate. However I am unable to successfully map a image created from QGIS to the sphere. The map appears slanted to one side.

I assumed this was because the spherical textures are different from the images produced by QGIS. Now using QGIS and https://www.naturalearthdata.com/ for maps would be a great asset for my project.

Is there a way to map this into the texture format that fits my OpenGL sphere?

Update: The map created in QGIS can be applied the Miller Cylindrical Projection. This allows mapping of the natural earth data (which is quite good) to an OpenGL texture map (according to my code). However this map becomes increasingly inaccurate as the poles are approached with respect to latitude.

If any one has any other suggestions, I'm still working on a better solution.


texts: Polygonizer Error "'NoneType' object has no attribute 'isMultipart' See log for more details"

texts:
I'm trying to polygonize. But I'm getting an error message:

  'NoneType' object has no attribute 'isMultipart' See log for more
  details

Basically what I'm trying to do is divide the city into blocks so I can assign a value to each block.

I'm not sure if I'm getting the error because of something about the lines or someting in the attribute table... Or something else.

Here's a closeup of what the street centerlines look like.

Maybe I need to clean up the lines or something, not sure. But anyway here's the error message I get when I run it.


texts: ArcPy Dissolve Statistics Fields - fixed into variable

texts:
How can I change this fixed POP column into my variable: field = parameters[1].valueAsText. It says syntax error everytime.

StatField = "POP MIN;POP MAX"

arcpy.Dissolve_management(in_features, region_Dissolve, "Class", StatField, "SINGLE_PART", "DISSOLVE_LINES")


texts: Attempt to generate ESRI REST Token generate 'no client_id specified' error even when client_id is in the request body

texts:
I'm making a HTTPS POST request to the following URL: http://www.arcgis.com/sharing/rest/oauth2/token/

In the Request body I've added the following JSON:

{"client_id":"<ID>","client_secret":"<Secret>","grant_type":"client_credentials","expiration":"1440","f":"json"}

Despite the client_id being clearly specified I always get the same response to my request:

{"error":{"code":400,"error":"invalid_request","error_description":"client_id not specified","message":"client_id not specified","details":[]}}

Is there something about the REST call I'm missing? Should I still be using JSON in the request body?


texts:
I'd like to incorporate the VIIRS Snow Cover Data Product (global
version VNP10C1) into my GEE code because of it's increased 375m-resolution compared to the MODIS Snow Cover Product (500-m resolution).
The only workflow I can find is to upload the VIIRS Snow Cover Data Product into ArcGIS for my dates of interest (2018-2020) and area of interest (Scandinavia), convert to raster (Band CGF_NDSI_Snow_Cover), then upload to GEE as multiple .tiffs.
Is that the best method for dealing with data for multiple years and an aoi of that size? Is there a way I can link to the data (https://nsidc.org/data/VNP10A1F/versions/1) from GEE and pull directly to GEE so I can visualize it there?


texts:
I need to create a query where if a property with the same lotplan has the same address it is only displayed once.

In case 1 we have the same lotplan and I want it to just show 32 once.

In case 2 -  I want just 63-65 to show on each lot

In case 3 it should just show 60 once.

The labeling query I have at the moment is using the following. Ideally the answer to this should incorporate these conditions as well.

Left$(ST_No_Label, Len(ST_No_Label)*Int(ST_No_Label<>"0"))+Left$(Chr$(13)+LOTPLAN, Len(LOTPLAN)*Int(Suf_No_1="R")*Int(ST_No_Label<>"0"))

See the following for background info If-then expression in Mapinfo label expressions

Without incorporating the above into a query we end up with a map that looks like this.

We can't just turn off duplicate labeling as it then doesn't label other lots with the same number but on a different street.

Can I do this without having to create a new column that has only the data I need? I would rather not do this as the backend DCDB changes often.


texts: Difference between ArcPy and ArcGIS API for Python?

texts: Could someone explain, from a high level overview, the difference between ArcPy and the new ArcGIS API package in Python?


texts: Transfering raster attributes to a shapefile

texts: I am using Arc GIS 10.2.  I have a raster layer and a shapefile layer. I would like to transfer the attributes of the raster layer to the shapefile layer using a majority rule. The raster layer is a lower resolution version of the shapefile layer with different attributes that I would like to have included with the higher resolution shapefile. I want to keep the attributes of the polygons in the shapefile layer and add the attributes from the raster shapes that make up the majority of each shapefile polygons.


texts: Insert cursor raises exception on geodatabase

texts:
I'm having trouble to insert a line geometry into a personal geodatabase (using Python/ArcGIS 9.3.1). Everything works fine if the target feature is a shapefile. My guess is that calling:

traj.shape = lineArray

is causing the exception:

 File "C:\Python27\lib\site-packages\win32com\client\dynamic.py", line 570,
 in __setattr__ raise AttributeError("Property '%s.%s' can not be set." %
 (self._username_, attr))AttributeError: Property '<unknown>.shape' can not be set.

However, lineArray holds the coordinates of the vertexes. To populate a new row I call:

newTraj = rowInserter.newRow()
newTraj.shape = lineArray
rowInserter.insertRow(newTraj)
lineArray.RemoveAll()
del rowInserter

Is there maybe someone with some experience on this topic?


texts: How to create Feature Class based on DBF schema?

texts: I have an attribute table schema (about 30 fields) stored in an Excel document, which I exported as a dbf. I can easily make a geodatabase table from this, using the Excel sheet as a template. How can I create a feature class from this schema (either from the Excel sheet or the dbf or the geodatabase table)? The Excel sheet is blank except for the field names.


texts:
Let's say I have the following geographical point represented in GeoJSON:

{"type":"Point","coordinates":[-109.0,35.292506]}

When I run the following conversion, things work as expected:

SELECT ST_AsText(ST_GeomFromGeojson('{"type":"Point","coordinates":[-109.0,35.292506]}'));
   st_astext
-----------------------
 POINT(-109 35.292506)
(1 row)

The problem occurs if the string passed to ST_GeomFromGeojson() omits the trailing decimal on -109.0:

SELECT ST_AsText(ST_GeomFromGeojson('{"type":"Point","coordinates":[-109,35.292506]}'));
   st_astext
----------------
 POINT(-109 35)
(1 row)

Notice how in this case rounding to an Integer value occurs on both latitude and longitude. The loss of precision on the latitude value in this case is very significant (32.53 km).

Note: I left out the CRS field for clarify, the problem occurs either way.

The same problem does not occur using the ST_GeomFromText() function. Both numbers are not rounded to Integers if the first looks like an integer.

My question is this: Am I missing something regarding the GeoJSON conversion function, or is this really unexpected behaviour? If it is a bug, where should I report it? Thanks!


texts: Determine if two geometries are similar in shape and location

texts: Suppose I have two geometries that have almost identical location and shape. I say almost on purpose here, since one or the other might have an extra node, or the order of coordinates in GeoJSON might differ. But is there a way to still determine their similarity, using either @turf utilities or PostGIS methods, perhaps within a certain error margin?


texts: Merge partial (or entire) geographical paths

texts:
I have data of a vehicle route and some data already inserted in the system. I have to find the parts of the route that are already presents in the system in order to avoid duplicate the info.
The data is given in latitude-longitude format. The positions given are never the same.
What do you think is the best method to make a comparison and determine which parts of the routes are equal enough?
Maybe the next image helps:

In yellow the parts of the rowute that have been marked as "similar"


texts: postgis - wrong ST_Distance with geography points

texts:
I'm trying to calculate the distance between to geographical points, but the calculated distance is bigger than the real one (from Google Maps).

The location column in postgresql is define as geography(Geometry,4326), and the query is:

  select ST_Distance(u1.location, u2.location)

For example, I have a pair of points for which the distance is:

160 km on Google Maps
190 km calculated with the above query.

What can be wrong? I may be missing a projection matter.

I have also tried the following query, but I get the same results:

  select ST_Distance(ST_GeogFromText(ST_AsText(u1.location)),
  ST_GeogFromText(ST_AsText(u2.location)))


texts:
I am trying to do a simple calculation that populates a field based on a set of acreage figures.  I am using if/else if statements with acreage ranges.  Here is the python script I came up with.

def reLabel(acres):
  if acres > 0 and acres <= 10:
    return "1"
  elif acres > 0 and acres < 20:
    return "2"
  elif acres> 0 and acres < 30:
    return "3"
  elif acres > 0 and acres < 40:
    return "4"
  elif acres > 0 and acres < 50:
    return "5"
  elif acres > 0 and acres > 50:
    return "6"
  else:
    return "0"

relabel (GIS_acre)

I hope that shows up correctly...

I keep getting an error that says syntax error unexpected EOF while parsing (<expresion>, line 1)

Any ideas where I'm going wrong?


texts: Leader lines not rendering when converted to PDF

texts:
This is very puzzling, I have an annotation feature class stored in a geodatabase that contains labels with leader lines.

The leader lines show up perfectly in the .mxd, it is only when I export to PDF that there is intermittent rendering of the leader lines. The puzzling part is that it works for some leader lines and not for others. The actually text of the annotation shows up in every single instance. A portion of the leader line stemming from the text shows up in every instance as well, it is only whether or not the leader line is completed that seems to be the symptom.

Please let me know if you need anymore info to help diagnose the problem.

Here is an example of the leader line on the left not being completed and the leader line on the right rendering correctly:


texts: Printing OpenLayers map with vector tile layers has CORS error for GeoServer vector tiles

texts:
I want to create a client side printing module for OpenLayers (v6) map. In OpenLayers web site there is a printing example with JSPDF (Export PDF example
). I used this example but there is an error:

In firefox:

SecurityError: The operation is insecure.

In Chrome:

SecurityError: Failed to execute 'toDataURL' on 'HTMLCanvasElement': Tainted canvases may not be exported.

and error is on this line:

mapCanvas.toDataURL('image/jpeg'), 'JPEG', 0, 0, dim[0], dim[1])

My Problem is that when I add Vector Tile Layers (not WMS, not vector layer, not OSM Tiles)  with MVT() format, layers are shown but after using print module, this error is showing in console.

I used crossOrigin: 'anonymous' like this answer for my vector tile source but the mentioned error there is on console and it is not working truly:

new VectorTileSource({
        crossOrigin: "Anonymous",
        format: new MVT({
            idProperty: "id"
         }),
        renderMode:"hybrid",
        defaultDataProjection: "EPSG:900913",
        url: GEOSERVER_ADDDRESS+"gwc/service/tms/1.0.0/" + STORE+":"+LAYER_NAME +
                "@EPSG%3A"+"900913"+"@pbf/{z}/{x}/{-y}.pbf"
})

My GeoServer is running on port 8080 and is not installed on Tomcat. Is there any solution for this issue on client side (OpenLayers) or server side (GeoServer)?

I couldn't find any solution for my case and I don't know how should I probably config GeoServer to solve this problem.


texts: QGIS 2.6 SAGA Processing Algorithm, "Add grid values to shapes", produces NULL Values

texts:
I'm using QGIS 2.6 Brighton on Windows with SAGA 2.0.8 installed using OSGeo4W.  I have the same versions installed on OSX as well which is not having this problem.

On Windows, I'm trying to extract elevation values from a DEM using a polyline shapefile by manually running the SAGA "Add Grid Values to Shapes" processing algorithm.  I've done this operation many times with the same data set on QGIS 2.4.  The DEM and shapefile are available as a zip file here: https://code.google.com/p/geospatialpython/source/browse/road.zip

The DEM and the shapefile are the same projection (EPSG:26910).  Everything runs without errors and the result is loaded, but what would normally be the averaged elevation value in a newly generated attribute column is NULL.  I can do the same thing on my OSX installation and get the expected result.  My end goal is actually a much more complex Python script that was working on QGIS 2.4 so this example is a trivial troubleshooting step to figure out where the problem is.

The algorithm command in the log is the following:

processing.runalg("saga:addgridvaluestoshapes","C:/qgis_data/road/road.shp","C:/qgis_data/road/dem.asc",0,None)

And there are no warnings or errors either in the Processing Toolbox log or in the log window of the tool's dialog.

Any ideas?


texts: Incorrect date displayed with Esri JavaScript API v.2.6

texts:
I am working on a site that is using version 2.6 of the Esri Javascript API.  I have an info window that displays some date information, and the date is one day behind the date in the database.  For example:

If I identify a feature in ArcMap, the date returned is 8/15/2007.
When I query for the same feature through the Rest endpoint I get 2007/08/15 00:00:00 UTC.
When I identify the feature in the web application, the date returned is 8/14/2007.

In the javascript code I am formatting the date as follows:

${SALEDT:DateFormat(selector: 'date', fullYear: true)}

What is causing this, and what can I do to have the correct date returned?


texts: Simplify two separate shapefiles, but ensure common edges are identical

texts:
I have two shapefiles.

One shows the political regions of Scotland.

The other shows the political constituencies of Scotland.

I need to simplify them to display in a web page. The problem is that when I simplify them there are small discrepancies between edges on the regions compared with the same edges in the constituencies.

Is there any known tool or technique I could use to allow me to do the following:

Simplify the regions shapefile and the constituencies shapefiles
Ensure all edges which are common to both are the exact same
Export a simplified regions shapefile and simplified constituencies shapefile

Up until now I've been trying it with Mapshaper.org and experimented with the different simplification algorithms and ratios. But I've not found a way to make it sync common edges between both. I'm aware of QGIS and ArcGIS, but have no idea if this kind of thing would be possible in them.


texts: Creating a grid from pixels of a georaster in Python with GeoPandas

texts:
Let say we have a georaster file with a given resolution. Is is possible to create a grid where each cell of the grid correspond to each pixel? Is it possible to store this grid as a GeoPandas dataframe?
This is a example of the image that I have.
fname = 'myFile.tif'
f = rasterio.open(fname)
show(f)

This is what I tried to do but I am not able to convert the meshgrid as a geopandas dataframe.
from affine import Affine
import xarray as xr
da = xr.open_rasterio(fname)
transform = Affine.from_gdal(*da.attrs['transform'])
nx, ny = da.sizes['x'], da.sizes['y']
x, y = np.meshgrid(np.arange(nx)+0.5, np.arange(ny)+0.5) * transform

The goal is to have a GeoPandas dataframe in which each rows correspond to a pixel with the value of the georaster because I want to make spatial join with other GeoPandas dataframes.


texts: "Cannot find variable: staticmethod" message while running Processing script

texts:
Every time I run my own processing script in QGIS 2.18.22 I get this message fifteen times in Processing tab in Log Messages Panel :

  1  Cannot find variable: staticmethod
  1  Cannot find variable: staticmethod
  ...

I don't get any error message and the script works fine (at least I think so). After a while I noticed that my processing script has fifteen staticmethods I defined in different classes. Numbers of the message increase when I add another staticmethod.

The outline of the script:

# Get parameters

# Import related packages

def RUN():
    pass

class TOOL():
    @staticmethod
    def Remove_Layers():
        pass

    @staticmethod
    def Clear_Map_Canvas():
        pass

class PROCESS():
    @staticmethod
    def Centroids(points, field):
        pass

    @staticmethod
    def Point_Grid_In_Region(gb, nat):
        pass

RUN()

Why do I get this message and could it be affecting results of the script?


texts: How I can use the pop-up window in OpenLayers on mobile? I have implemented the FramedCloud popup on my website, but when I open my web browser on a mobile device, the popup window does not seem to work or listen to any open action open.


texts: Installing Apache Web Server from OSGeo4W

texts:
I am trying to create a QGIS web client on my Windows 8 But with OSGEO4W not appear the option to make Apache installs Web Server, as of this writing in several tutorials and documentations.
A print screen of my Windows 8

Should this well, this print belongs to a tutorial found on the web


texts:
I have a layer in GeoServer and its projection is EPSG:27700. When I query the layer in my app I get and error.

Explanation:

The GeoServer preview displays it correctly. So I copy the JavaScript from the OpenLayers preview page of the layer and use it in my map, and then I compare the browser console for network requests. My map sends an SRS parameter of EPSG:3857 in its requests, whereas the GeoServer preview sends an SRS parameter of EPSG:27700 in it's requests. When I copy the request and past it into the browser I get a WMS file with this seemingly unrelated error message:

Error rendering coverage on the fast path
java.lang.IllegalArgumentException: The reference geometry must be either Polygon, MultiPolygon, or a polygonal PreparedGeometry
The reference geometry must be either Polygon, MultiPolygon, or a polygonal PreparedGeometry

If I manually modify the SRS value in the URL and paste it into the browser again... I get a valid png file.

I assume the SRS it is sending is identifying what the projection of the map object and I cannot change it? Can anyone help me to understand what is going on please?

=== marked as duplicate, although the question linked doesn't help me much ===
By "copy the request" I meant that Using the browser console I copied the URL my map was firing to GeoServer to request tiles.

The only difference between what my map's URL and the URL from the GeoServer's OpenLayers preview is that the SRS parameter is different.

The GeoServer preview request gets a valid png, the request from my map results in the error mentioned. I guess specifically I want to understand how that results in the error as I cannot understand how referencing a geometry even comes into this scenario.


texts:
I want to know what is the official survey datum, and projection system for Zimbabwe? I think it is ARC1950, UTM Zone36.

Is it better to convert all the old coordinates to WGS LO29 (LO depending on nearest uneven longitude)?


texts: Delimited text layer doesn't align with shapefile

texts:
I have a shapefile of 25 major drainage basins of Canada and a .csv file of 77 point data I would like to plot within these drainage basins. I am using QGIS 3.10.9
I reprojected the drainage basins from WGS83 to NAD83/Statistics Canada Lambert EPSG:3347 because I need the measuring unit to be in meters so I can create buffers around point data. The point data are saved as a .csv and I am trying to add them as a delimited text layer (and afterwards will export to shapefile).
When I specify the same CRS (EPSG:3347) and then click ADD, the points show up in a completely different area (and seem way too close together).
I have tried to save as a shapefile, reproject, define projection. None of these seem to work. The points remain off to the left.
Here are some screenshots so you can see - anyone know where I might have gone wrong?
Here is the point data: https://www.canada.ca/en/environment-climate-change/services/environmental-indicators/water-quality-canadian-rivers.html


texts: QGIS Expression - Show multiple records of other map layers

texts:
I am working of a map layout in QGIS 3.10. In my layout i am inserting a HTML block. In that HTML block i want to show content by building an expression. This expression will show information from 3 different layers present in map. In below image as show Fields and Values is showing attributes of one layer which is being used in map layout Atlas. I want to utilize other layer attributes as well. THere is an option to to do that so i used following

 get_feature( 'division1_8cc9f2a7_078a_4d49_8e5a_826e54bc33f1', 'memberOr_5',   "STATE_ABBR"  )

But this is returning me first feature of layer which satisfies the queries but there are more then one feature how can i get others ?


texts: ReduceRegions Graph JavaScript in Google Earth Engine

texts:
I'm trying to get NDVI values as a featureCollection with my nodejs server. But the "ReduceRegions" function does not recognize the input. How to solve this problem?

EE code link: (where it works fine) https://code.earthengine.google.com/3f8355e66a11f2a4db69d82b350deb90

  Error: Unrecognized argument type to convert to a FeatureCollection: [object Object]
      at new ee.FeatureCollection (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:16907:17)
      at Function.ee.promote_ [as promoter_] (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:17822:53)
      at ee.ApiFunction.ee.Function.promoteArgs (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:15799:40)
      at ee.ApiFunction.ee.Function.apply (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:15776:49)
      at ee.ApiFunction.ee.Function.callOrApply (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:15793:15)
      at ee.Image.destination.(anonymous function).signature.destination.(anonymous function) [as reduceRegions] (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:15945:24)
      at reduce (C:\Users\Omen-\Desktop\heroku\ams-ocp\functions\graphNdvi.js:25:22)
      at C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:16505:12
      at new ee.CustomFunction (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:17444:24)
      at Function.ee.CustomFunction.create (C:\Users\Omen-\Desktop\heroku\ams-ocp\node_modules@google\earthengine\build\main.js:17502:10)

var app = function (image) {
        var ndvi = image.normalizedDifference(['B5', 'B4']);
        image = image.addBands(ndvi.rename('NDVI'));
        return image;
    };

var collection = ee.ImageCollection('COPERNICUS/S2')
        .filterBounds(inputs.reg)
        .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
        .filterDate(inputs.date[0], inputs.date[1])
        .map(app);

var reduce = function (image) {
        return image.reduceRegions(
            {
                collection: inputs.reg,
                scale: 30,
                reducer: ee.Reducer.mean()
            })
}

var ndvi = ee.ImageCollection(collection.select(["NDVI"]))
var featureCollection = ndvi.map(reduce).flatten();


texts: Convert TIFF image (Mosaic Dataset) to Polygon Feature

texts:
I have geo-referenced TIFF images, I need to convert them in polygon feature class. I have used Raster to Polygon Toolbox from ArcGIS but it ran for more then 48 hours and did not give me meaningful data.

I need some suggestion on this so that I can get all polygon features from TIFF images.  My images are scanned copy of polyvinyl sheets.


texts: How to determine if a polygon is self-intersecting using either ArcPad 8.0 or ArcMap 9.3.1?

texts:
I am try to automate a sync replica mxd a nightly basis.  The issue I am running into is that in some the of the AXFs (that I get from the users and check-in) have polygons that are self-intersecting and is causing problems with the sync replica.  The sync errors out and the polygons need to be fixed manually.  Is there any way to detect if a polygon is self-intersecting in ArcPad 8.0, using VBA or some other technique, to stop the creation of the polygon?  If not, is there some tool in ArcMap 9.3.1 that can be run again the mxd to detemine the self-intersecting polygons and fix them?  Any ideas will be greatly appreciated.

Thanks.

Tim


texts: How to use gdal utilities to subset from a raster?

texts:
How to you use gdal to subset a smaller raster from a bigger one given the bounding coordinates?

EDIT: Can you subset a region from a vrt ?


texts:
Does Anyone have legend colors for the NHD Flowline that USGS sets colors for theirs ?

I'd like to get a copy of their color legend.

They can be from their website, brouchure, lyr files or anything that I can set it up easily to match their colors.


texts: InsertCursor in ArcGIS Pro doesn't understand my coordinates

texts:
I am trying to use the arcpy.da.InsertCursor to add simple points to a feature class, but it doesn't work for me as every single point I add ends up somewhere near Paris. Here's the simplest code block that doesn't work for me:

arcpy.CreateFeatureclass_management(' (my path here) ...\Map.gdb','Hotspots2', 'POINT')
a = arcpy.Point(7, 7)
cursor = arcpy.da.InsertCursor("Hotspots2",["SHAPE@XY"])
cursor.insertRow([a])
del cursor

It ends up here:

at coordinate ~ (3.3, 48). Whatever coordinates I enter into a don't matter for the outcome. If I change a to a point geometry it doesn't work. If I use SHAPE@ and a point geometry it also doesn't work. Every darn point ends up there. What am I doing wrong and can someone test to see if they get the same result??


texts: Using HERE API credentials for Batch Geocoding?

texts:
I have been trying to create my own geocoder using the HERE API to handle geocoding jobs. I have the HERE api I'D, code, and key. I have reviewed much of their developer information and am still confused about how to use the developer material to geocode batches of addresses. Some material appears as though I click on a url (such as https://batch.geocoder.api.here.com/6.2/jobs) while the documentation seems show chunks of script.

How do I use my HERE API credentials to actually geocode multiple addresses?

Do I need to create a website using HTML so that I can create a place or web portal to upload my data as in the example of Batchgeo or other web-based geocoders?

That was what I was attempting initially.

How do I use the HERE developer examples and materials to actually geocode a list of multiple street addresses in the U.S. and receive either a geocoded file (.shp or GeoJSON) that I can use in a map?


texts: Inset map in print composer shifts when saved as pdf

texts: I'm working on a map in the print composer window of QGIS, with an inset map showing a larger extent and highlighting the area of my main map as an overview. Everything looks good in the print composer, but when I save as a pdf or image, or print as a pdf, the inset map stays the same in the print composer but the inset map on the pdf has totally shifted position and scale. I have the layers and styles for the map item locked.


texts: D maps using Openstreetmap data with the help of Qt. I need to use Qt specifically as other components in the software have objects that integrate into the map using Qt. How should I proceed? Can I parse postgis data to construct maps using Qt? Are there some tutorials?


texts: Using VBscript restrict layer editing in ArcPad

texts:
I am attempting to restrict editing on certain layers based on the user log in.

Dim objLyr1
Dim objLyr2
Dim objLyr3
Set objLyr1 = Map.Layers("PatchPanel")
Set objLyr2 = Map.Layers("wHydrant")
Set objLyr3 = Map.Layers("ElectricPoles")
If UID = 2 Then
 Msgbox "Electric user. Only enable editing for Electric Users."
 objLyr1.Editable = True
 objLyr2.Editable = False
 objLyr3.Editable = True
ElseIf UID = 4 Then
 Msgbox "Water user. Only enable editing for Water Users."
 objLyr1.Editable = False
 objLyr2.Editable = True
 objLyr3.Editable = False

ElseIf UID = 3 Then
 Msgbox "Sewer user. Only enable editing for Sewer Users."
 objLyr1.Editable = False
 objLyr2.Editable = False
 objLyr3.Editable = False

Else
 Msgbox "Invalid user. Editing Unavailable"
 objLyr1.Editable = False
 objLyr2.Editable = False
 objLyr3.Editable = False
End If

This all happens after the log in is verified and the if statement works. The error message I'm getting says "error object required: 'objLyr1' Source Text Unavailable" It is error code 800A01A8.

I'm thinking the Map.Layers command may be generic for instructional documentation and I need to tailor it to my program but that's just a guess.

I'm not sure why and it is a rather vague error.

Does anyone have a clue to why I'm receiving this error?


texts: Where can I find shapefiles or geospatial data for postal codes/ Código de Endereçamento Postal for Brazil

texts: I am looking for shape files for Brazilian postal codes. I have looked at IBGE site and other sources but could not find anything. Free shape files will be better since we plan to use them for commercial purposes.


texts:
I have two layers, one of points (locations of fruit trees) and one of Atlanta zip codes, in ArcMap 10, but the layers are not aligning. The fruit trees are in a corner and the map is in the middle top.

The fruit points were successfully converted to a shp file from a csv file.

The two layers have the same projected coordinate system NAD_1983_StatePlane_Georgia_West_FIPS_1002_Feet

and the same geographic coordinate system
GCS_North_American_1983

I thought the error would be that the systems didn't match up, but that's not the case.


texts: How can I select the Text Annotation tool through python

texts:
Some of the QGIS tools can be selected/activated via PyQGIS using the iface object (see QgisInterface). For example:

iface.actionIdentify().trigger() # Activates the Identify tool

However, the Text Annotation tool (and also HTML, SVG, and Form annotation tools) are not available in iface.

How can I select/activate the Text Annotation tool from my QGIS python plugin?


texts: Chrome not showing images in pop ups in ArcGIS Online

texts:
Google Chrome is not showing an image in an ArcGIS Online pop-up.
Is there any way to fix this problem?


texts: Difficulty with using truetype font symbols in Mapserver

texts:
I am basically following this workflow and having no luck. I have add a truetype font file to my fonts.list (under the alias 'ddv') however, I cannot get the fonts to render in the service. Any ideas, maybe...?

    # ma_interstate_labels.map

########## English ##########
  LAYER
    NAME "county_labels"
    GROUP 'county_labels'
    TYPE ANNOTATION
    STATUS DEFAULT
    DATA "ROADS_PT_MA_Interstate_3785.shp"
    LABELITEM "RTE_NUM1"
    # Visible in map from zoom level 6 onwards
    MAXSCALEDENOM 10400256
    MINSCALEDENOM 2539
    # shell symbol
    SYMBOL
      NAME shell
      TYPE truetype
      FONT "ddv"
      CHARACTER "&#35;"
    END
    CLASS
       STYLE
         SYMBOL shell
         ANTIALIAS TRUE
       END
      LABEL
       TYPE truetype
       FONT "deja"
       SIZE 5
       COLOR 52 88 133
       ALIGN center
       WRAP " "
       MAXLENGTH 4
       ANGLE 0
       PARTIALS FALSE
       END # end label
    END # end class
  END # end of layer


texts: I have a problem with saving as GMT in QGIS. There is no option in my QGIS version to save as in GMT. How I can add this option in QGIS?


texts: Add Cartoview apps to an existing GeoNode installation

texts:
I´m using GeoNode 2.6.3, running on ubuntu server 16.04. I´ve seen some Cartoview apps that might be usefull for our porpouse, but i´m not sure if this instructions work on a GeoNode instance that is alredy running. It instructs to install (again) postgresql, geoserver, and other that GeoNode is alredy using.

Will it work if i skip the first steps and configure the new databases on the existing Postgresql 9.5?
If I install cartoview, is it safe for all exsisting data (layers, maps, users, groups) to keep working fine?

If any has done so, please confirm that it works well. I just don´t want to start the instalation and mess up what is alredy working fine.


texts: QGIS/QDAL: Accessing OGR_STYLE across all layers in a KML file

texts:
I get a regularly updated KML file which I need to pull into QGIS.
It has different folders/layers (about 12), and features in each folder/layer have different <styleURL>s. I need to access each feature's styleURL (along with name, description, and the location) for processing in QGIS.

1st attempt - just import the KML file inside QGIS. I get one layer per folder (well, lines and points separately), but the styleURL is not an accessible field.

2nd attempt by manually running ogr2ogr and using the special field OGR_STYLE:

ogr2ogr -f GPKG -oo LIBKML_RESOLVE_STYLE=YES -sql "select *, OGR_STYLE from Boundaries" twinrock_test.gpkg twinrock13.kml

Works fine for a single named folder, Boundaries, but I need to run it for all folders in the kml file, as if I just did ogr2ogr -f GPKG twinrock_test.gpkg twinrock13.kml (without the mucking around with OGR_STYLE).

3rd attempt:

ogr2ogr -f GPKG -oo LIBKML_RESOLVE_STYLE=YES -select Name,description,OGR_STYLE twinrock_test2.gpkg twinrock13.kml

Fails since the special field OGR_STYLE seems only accessible via -sql, not -select.

Any ideas?


texts:
I have a python script that works (it's worked for several other projects in the past and it worked for my current project initially). But whenever I run it, Python crashes sporadically. It would go through some iterations (I'm calculating zonal stats as table for 540 files for 10 countries, so over 5000 files in all, it's worked fine for the first country and now it's not working).

I'm not getting an error message or anything, it's just stopping. I mean I would get the:

      ========= RESTART =================

but that's it. Then I would click back to the script, click F5 and restart it, it would stop randomly and I'll restart it again. I have no idea why it's doing this and there's no way that I can babysit the program.

HELP?!


texts: Calculate mean values for all bands and date to a dataframe

texts:
I have 5 rasters, each having 10 bands, where each band corresponds to a specific date.

How do I calculate mean values for all 10 bands and save them?

I need to save the mean values in a dataframe like this example:

   rasters   bands       mean         date             years
    raster1   band1       1234         24/01/2000       2000
              band2       45213        01/03/2000       2000
              band3        4221        12/04/2000       2000
              band4         ...         ...              ....
              band10
   raster 2   band 1       4521        01/02/2001      2001
              band2        1234        04/05/2001      2001
              ...           ....         ....

This script calculates mean value for each band in my raster:

raster = gdal.Open("D:/script/NDVI2000.tif")
bands = raster.RasterCount

for band in range(1, bands+1):
    data = raster.GetRasterBand(band).ReadAsArray().astype('float')
    mean = np.mean(data[data != 0]) #calculate mean without value 0
    print("Band %s: Mean = %s" % (band, round(mean, 2)))
Out[2]:
Band 1: Mean = 1712.83
Band 2: Mean = 1803.14
Band 3: Mean = 1662.33
Band 4: Mean = 1868.77
Band 5: Mean = 1900.97
Band 6: Mean = 2031.13
Band 7: Mean = 1847.89
Band 8: Mean = 2185.66
Band 9: Mean = 2081.14
Band 10: Mean = 8248.7

How can I save the mean result stats in the dataframe using gdal python 3 after I have the results for all rasters?

I need this dataframe to plot times series for NDVI.


texts: How do I Calculate Grouped Percentiles or Batch percentiles in ArcMap?

texts:
My question is similar to this one which was answered successfully: Calculating Percentiles in ArcMap?

My addition is how do I batch this solution for multiple files or apply a group-by feature to this code?

I have a large dataset with HUC14 watershed classifications and groundwater recharge values. I would like to determine the quintile that each record is in for each HUC14. So, first the data needs to be grouped by HUC14 and then the percentile calculated. Right now I have the data as a single feature class (feature layer) and also as separate files (shapefiles), so either a batch solution or a grouping solution would be ok.

EDIT:
I wanted to add the code that I finally used based heavily on @Farid Cher 's response.

import arcpy
import numpy as np
import os

#loop through all Shapefile in a folder and call the CalcPercentile method

workspace = "X:\Jocelyn\GWR\HUC14b"
walk = arcpy.da.Walk(workspace, datatype="FeatureClass")
for dirpath, dirnames, filenames in walk:
    for filename in filenames:
        featureClass = os.path.join(dirpath, filename)

    #First add the percentile Rank Field
         arcpy.AddField_management(featureClass, "PercRank", "LONG", 5, "", "","", "NULLABLE")

        inputFeatureClass = featureClass

    #Creates a feature layer to allow for the selection and removal (using switch selection) of records that have a GSR code of 999 (water and wetlands) or no soil group associated with it.
    #These records should not be used when calculating the percentile rank.
        FeatureLayer = arcpy.MakeFeatureLayer_management (inputFeatureClass, "temp_lyr")
        FL_Selection1 = arcpy.SelectLayerByAttribute_management (FeatureLayer,"NEW_SELECTION", """ "SoilGroup" = ' ' OR "GSR32_code" = 999""")
        FL_Selection2 = arcpy.SelectLayerByAttribute_management (FL_Selection1,"SWITCH_SELECTION")

    #Only uses the selected features in the CalcPercentile function
        CalcPercentile(inputFeatureClass)

    #Deletes the temporary feature layer to avoid cluttering and space issues.
        arcpy.Delete_management(FeatureLayer)

def CalcPercentile(inputFeatureClass):
    arrp = arcpy.da.FeatureClassToNumPyArray(FL_Selection2, ('GWR_V'))

    arr = np.array(arrp,np.float)

#to create 5 ranks
    p1 = np.percentile(arr, 20)  # rank = 1 (lands that provide the lowest volume recharge within the HUC14)
    p2 = np.percentile(arr, 40)  # rank = 2
    p3 = np.percentile(arr, 60)  # rank = 3
    p4 = np.percentile(arr, 80)  # rank = 4
    p5 = np.percentile(arr, 100)+1  # rank = 5 (lands that provide the highest volume recharge within the HUC14)
#Print the quintile breaks that are calculated.
    print "p1=%s" % p1
    print "p2=%s" % p2
    print "p3=%s" % p3
    print "p4=%s" % p4
    print "p5=%s" % p5

#use cursor to update the new rank field
    with arcpy.da.UpdateCursor(inputFeatureClass , ['GWR_V','PercRank']) as cursor:
        for row in cursor:
            if row[0] < p1:
                row[1] = 1  #rank 1
            elif p1 <= row[0] and row[0] < p2:
                 row[1] = 2
            elif p2 <= row[0] and row[0] < p3:
                 row[1] = 3
            elif p3 <= row[0] and row[0] < p4:
                 row[1] = 4
            else:
                 row[1] = 5

        cursor.updateRow(row)


texts: How to solve a network error with a WMS?

texts:
There is an error with a WebMapService I am using. The error message is depicted in the following image:

The properties of the file are shown in the image below.

Does anyone know how to fix this corrupted file? I am using QGIS version 3.4


texts: Can't select vector layer QGIS Zonal Statistics

texts: I'm running QGIS 3.4.0 on Ubuntu 18.04.1.  I have a raster layer and a vector layer from a .shp file in the same CRS.  When The window pops up for zonal statistics, I can select my raster layer but when I click the dropdown box to select my vector layer, nothing happens.


texts: Efficiency of reading non-tilled geotiffs in GDAL?

texts:
I need to stream large geotiff files through web server.
By large I mean 1Gb-3Gb.

The datasets I was given aren not tilled (pixel strips).

I wonder how does it affect efficiency of reading them using GDAL.
If I request only part of image - does GDAL load into memory only the relevant part of image?
Does the file need to be tilled in order to be read piece by piece?

This is how gdalinfo output look like for typical file I want to read:

Band 1 Block=43554x1 Type=Byte, ColorInterp=Red
  NoData Value=-10000
Band 2 Block=43554x1 Type=Byte, ColorInterp=Green
  NoData Value=-10000
Band 3 Block=43554x1 Type=Byte, ColorInterp=Blue
  NoData Value=-10000
Band 4 Block=43554x1 Type=Byte, ColorInterp=Alpha
  NoData Value=-10000

Question is:
if i want to access 512x512 fragemnt will GDAL read (best case scenario)
at least 43554x512 pixels? (and allocate at leat that much memory).


texts: Sampling points in R with both distance constraint and group

texts:
I have a spatial point data frame with categorical variable as a column named "type", I would like to draw a random sample from these points (let's say 100 points) with maximum and minimum distance constraints but also grouped by type (i.e 20 points for each type supposing "type" has 5 unique values)
my data
class(my_data)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

I can get a sample with distance constraints using spatialEco package like this
library(spatialEco)
sample <- subsample.distance(my_data, size = 100, d = 0.5, d.max = 3)

or grouped sample using sample_n from dplyr after converting to data.frame like this:
library(dplyr)
my_data_df <- as.data.frame(my_data)
sample <- my_data %>%
                      group_by(type) %>%
                      sample_n(20)

I'm looking for a way that do both at the same time!


texts: How to install Sourcepole zonal statistics

texts:
I am trying to install the zonal statistics tool created/supplied by sourcepole from http://www.sourcepole.com/2010/12/16/zonal-statistics-plugin-for-qgis.

The instructions area as follows:

**To install: unpack zonal_statistics.tgz under /src/plugins, insert zonal_statistics as new subdirectory in /src/plugins/CMakeLists.txt and recompile.**

I have downloaded and unzipped the file, but I am now not sure how to add it to the QGIS software.  I am using a Mac OS, and QGIS version 1.7.4.

I can't find the folder /src/plugins which the instructions speak of.... and in general am very confused by the instructions (more due to my semi-computer illiteracy than the instructions themselves).

Any help I would greatly appreciate.

Regards

Becky


texts: How do you turn on the dark mode for QGIS 3.6.1? I see where people have an option in the general settings of the Project Properties menu, but I am not seeing that option in my version.


texts: Getting origin as long, lat in GDAL without knowing EPSG

texts:
My goal is to get the long, lat coordinates for the origin of a particular geotif.

gdalinfo outputs:
Upper Left (-2764486.928, 3232110.510) (130d 7'22.42"W, 49d56'48.81"N)
But I don't understand the meaning of (-2764486.928, 3232110.510) - for instance what units are being represented here?

What I want is something like:
Upper Left  (-125.0208333,  49.9375000) (125d 1'15.00"W, 49d56'15.00"N)
which gdalinfo yields on a different geotiff file that I have.

I'm actually using ds.GetGeoTransform() in python to get this same information gdalinfo outputs and am looking for a solution in python to convert (-2764486.928, 3232110.510) into a long, lat. I tried fmark's answer here, but it didn't convert (-2764486.928, 3232110.510) into a long,lat.

Here is the entire gdalinfo output

  Driver: GTiff/GeoTIFF
  Files: mask.tif
  Size is 2145, 1377
  Coordinate System is:
  PROJCS["unnamed",
      GEOGCS["Coordinate System imported from GRIB file",
          DATUM["unknown",
              SPHEROID["Sphere",6371229,0]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433]],
      PROJECTION["Lambert_Conformal_Conic_2SP"],
      PARAMETER["standard_parallel_1",25],
      PARAMETER["standard_parallel_2",25],
      PARAMETER["latitude_of_origin",25],
      PARAMETER["central_meridian",265],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]]]
  Origin = (-2764486.928100540300000,3232110.510093217700000)
  Pixel Size = (2539.703000000000000,-2539.703000000000000)
  Metadata:
    AREA_OR_POINT=Area
  Image Structure Metadata:
    INTERLEAVE=BAND
  Corner Coordinates:
  Upper Left  (-2764486.928, 3232110.510) (130d 7'22.42"W, 49d56'48.81"N)
  Lower Left  (-2764486.928, -265060.521) (121d33'48.72"W, 20d10'43.04"N)
  Upper Right ( 2683176.007, 3232110.510) ( 60d51'59.22"W, 50d 6'46.25"N)
  Lower Right ( 2683176.007, -265060.521) ( 69d11'55.57"W, 20d19' 6.55"N)
  Center      (  -40655.461, 1483524.995) ( 95d27' 9.15"W, 38d13' 5.62"N)
  Band 1 Block=2145x1 Type=Float32, ColorInterp=Gray

Here's what I tried in python:

geo_transform = ds.GetGeoTransform()
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())
wkt = """
    GEOGCS["NAD83",
    DATUM["North_American_Datum_1983",
        SPHEROID["GRS 1980",6378137,298.2572221010002,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6269"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4269"]]"""
new_cs = osr.SpatialReference()
new_cs.ImportFromWkt(wkt)
transform = osr.CoordinateTransformation(old_cs,new_cs)
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()

#get the coordinates in lat long
latlong = transform.TransformPoint(gt[0], gt[3])

The output is latlong = (-2764486.928, 3232110.510, 0.0)

Using BradHards's link, I've got the code pasted below and believe I just need to figure out the source geotiffs EPSG. gdalinfo says it's 9001, but when I try using that the code blows up with

TypeError: in method 'Geometry_Transform', argument 2 of type 'OSRCoordinateTransformationShadow *'

The geotiff originated as a GRIB file from nomads

ds = gdal.Open("mygeotiff")
gt = ds.GetGeoTransform()
long_lat = get_point_as_long_lat(gt[0], gt[3], ???what_epsg???)

def get_point_as_long_lat(coord_x, coord_y, input_epsg):
  # Spatial Reference System
  output_epsg = 4326

  # create a geometry from coordinates
  point = ogr.Geometry(ogr.wkbPoint)
  point.AddPoint(coord_x, coord_y)

  # create coordinate transformation
  in_spatial_ref = osr.SpatialReference()
  in_spatial_ref.ImportFromEPSG(input_epsg)

  out_spatial_ref = osr.SpatialReference()
  out_spatial_ref.ImportFromEPSG(output_epsg)

  coord_transform = osr.CoordinateTransformation(in_spatial_ref, out_spatial_ref)

  # transform point
  point.Transform(coord_transform)

  # return point in EPSG 4326
  return point.GetX(), point.GetY()


texts: Zoom hangs ArcGIS API for JavaScript application?

texts:
I'm using ESRI arcgis 3.13

I have map with minZoom set to 3.

When user tries to zoom in to level 19, the program just hangs and after a while, comes back with a generic browser error not from the application that says "more memory is needed" etc. So the application is killed during zoom. Shown here:

I tried to debug this, and pinpointed to this line, this is the line where the program would stop and hang, and it is calling an ESRI map method:

_map.centerAndZoom([-79, 43], 19);

I have this for my map intialization:

_map = new Map("myMap", {
    center: [-100, 48],
    zoom: 3,
    basemap: "streets",
    minZoom: 3,
    logo: false
});

However, this was all working before, and just stopped work this week.

Also, if I set initial map zoom to be 7, the app hanging problem goes away. So I am thinking if it because the zoom is changing too much. But like I said, the program was working before even with that much of zoom, and it worked before even from zoom level 1 to zoom level 19.

Any ideas how to fix this? Or where to look for potential problems?