texts: Extract features in .gdb file based on their attributes

texts:
I have a .gdb file with >17,000 features. Simplified example of attributes table:

ID      NAME    DATE    ETC
0001    london  2012    12
0002    london  2012    23
0003    thame   2012    45
0004    thame   2012    56
0005    kent    2012    67
0006    london  2011    78
...

I have an Excel file with a list of ~250 names:

NAME
kent
essex
sussex
london
...

I want to extract only the features with attributes that match the names in the Excel file and save them as shapefiles. Multiple features share the same attributes. How to do this using FOSS? I think there is a function in ArcGIS called FeatureClassToShapefile that does this.

I found this code to run GDAL directly from the terminal to convert the .gdb to a shapefile:

ogr2ogr -f "ESRI Shapefile" /test/extraction /test/test.gdb

I can select a subset of layers by appending these to the above code:

ogr2ogr -f "ESRI Shapefile" /test/extraction /test/test.gdb example_layer

If I try to use the same formatting to select a subset of features I get the error:

  FAILURE: Couldn't fetch requested layer 'example_feature'!

Is this possible using this method? I've looked at the GDAL documentation here but couldn't find an answer.

I'm open to other solutions, using GDAL, R pr QGIS preferably.


texts: GeoServer and NetCDF - Gap around prime meridian

texts:
I'm attempting to use GeoServer to publish NetCDF layers. So far has been working well using EPSG:4326 (the native srs) as reference system in the GetMap request. In this case the image from GeoServer looks something like below

(request parameters "bbox=-90%2C-77.5%2C90%2C79.5&width=768&height=768&srs=EPSG%3A4326&format=image/png")

If the reference system in the request is changed to EPSG:3857 (Mercator) however, a small gap appears around the prime meridian as seen in below image

(request parameters bbox=-10018754.17%2C-14106322.01%2C10018754.17%2C15225854.97&width=768&height=768&crs=EPSG%3A3857&format=image/png)

The NetCDF file has 1°x1° spatial resolution and longitude axis 0°-359°.

I have tried changing the request bbox parameters, the bounding boxes for the layer in GeoServer as well as changing the style to the default raster style, but the gap still shows.

How do I get rid of this gap?


texts: How do I detect when "Map data not yet available"?

texts:
My app uses the CartoDB static map API to get a "center" map at zoom level 19 or 20.  Most times it gets a good map, but depending where it's looking it sometimes gets a "Map data not yet available" map like this:

It uses the ESRI tiles http://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x} because these provide sufficient clarity for most of our coverage area (New Zealand).

I want to be able to detect when I get these (perhaps by running another command first to determine availability) so that I can use a different basemap instead, and avoid returning this kind of result to my customer.

This all happens within a C# program, where I can do all sorts of things programatically, but nothing manually.


texts: OSTN15 transform not working in postGIS build

texts:
I have an issue with my new (CentOS) postGIS setup not transforming WGS84 to BNG via OSTN15 correctly.  It all worked beautifully in postGIS2 with proj4.8. As recommended, the relevant entry is in the PostGIS spatial_ref_sys file, with a custom srid 27789 :
"+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +nadgrids=/usr/share/proj/OSTN15_NTv2_OSGBtoETRS.gsb"

ST_Transform(geom,4326,27789) and vice versa both worked as expected, picking up the gsb file..
However, With my latest PostGIS3 and proj7.2 build, the transform no longer works. It gives the standard transform rather than the OSTN15 one, as if it is no longer using the gsb file. gsb file permissions are identical between my builds, and the gsb file location has not changed.
Entering the custom proj string into the ST_Transform command directly works OK. For some reason it is not picked up from spatial ref sys anymore. Any ideas why please?
i.e. running the following code:
postGIS2 : "coords_expected" and "coords_via_OSTN15_spatialrefsys" are the same.
postGIS3 : "coords_expected" and "coords_via_OSTN15_spatialrefsys" are different -they ignore the OSTN15 transform.
postGIS3 : "coords_expected" and "coords_via_OSTN15_manual" are the same, so postGIS can read the gsb files.
--This to verify whether the OSTN15 transform is working or not...
    select uprn,
    x_coordinate::text || '_' || y_coordinate::text as coords_expected,
    round(ST_X(geomBNG_OSTN15_srs))::text    || '_' || round(ST_Y(geomBNG_OSTN15_srs))::text    as coords_via_OSTN15_spatialrefsys,
    round(ST_X(geomBNG_OSTN15_manual))::text || '_' || round(ST_Y(geomBNG_OSTN15_manual))::text as coords_via_OSTN15_manual
    from
    (   --do a OSTN15 transform on the latlon - one via spatial ref sys (srs) - other via 4326/OSTN15 params from spatial ref sys (manual)
        select t1.*,
        ST_Transform(geom4326, 27789) as geomBNG_OSTN15_srs,
        ST_Transform(geom4326, '+proj=longlat +datum=WGS84 +no_defs', '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +nadgrids=/usr/share/proj/OSTN15_NTv2_OSGBtoETRS.gsb') as geomBNG_OSTN15_manual
        from
        (
            --coords from OS ABP...  make latlon geom
            select '100090925668' as uprn,
            620537 as x_coordinate, 308016 as y_coordinate,
            52.6251053 as latitude, 1.2569831 as longitude,
            cast(ST_SetSRID(ST_MakePoint(1.2569831, 52.6251053),4326) as geometry(Point,4326)) as geom4326
        ) t1
    ) t2
    ;

    --output showing the OSTN15 spatialrefsys transform is being ignored
    "uprn"          "coords_expected"   "coords_via_ostn15_spatialrefsys"   "coords_via_ostn15_manual"
    "100090925668"      "620537_308016"     "620538_308019"             "620537_308016"

Other details :

--public.spatial_ref_sys entry (excluding srtext field for conciseness)
select srid, auth_name, auth_srid from public.spatial_ref_sys where srid in (4326, 27789);

"srid"  "auth_name" "auth_srid"     "proj4text"
4326    "EPSG"      4326            "+proj=longlat +datum=WGS84 +no_defs "
27789   "EPSG"      27700           "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +nadgrids=/usr/share/proj/OSTN15_NTv2_OSGBtoETRS.gsb"

--postgis build details
select postgis_full_version();

POSTGIS="3.0.3 6660953" [EXTENSION] PGSQL="120" GEOS="3.8.1-CAPI-1.13.3" PROJ="7.2.0" GDAL="GDAL 3.2.0, released 2020/10/26" LIBXML="2.9.1" LIBJSON="0.11" LIBPROTOBUF="1.0.2" WAGYU="0.4.3 (Internal)" RASTER

There are no errors reported, just an inaccurate answer (i.e. worst case scenario).


texts: Can I do the sql server for find fnode and tnode?

texts: Can I do the sql server for find fnode and tnode field?


texts: Displaying the cloud free composite of Sentinel2 using s2cloudless in GEE

texts:
I applied the tutorial  for 'Sentinel-2 Cloud Masking with s2cloudless' and have this script in GEE:
https://code.earthengine.google.com/c525b3926841d4953e016e64682a4e5c
In this tutorial, I do not understand the  purpose of 'cloud transform' layer and I do not know how to display the final cloud free image.
Should I apply a median to see the final result? If yes, over which layer?


texts: What are the minimum system requirements for QGIS 3.0 Desktop?

texts: What are the minimum system requirements for QGIS 3.0 Desktop, as speed of CPU and how many cores, how much disk space is needed for installation, size of HDD and swap space, minimum RAM memory, display and screen resolution, graphics memory, etc ?


texts:
I am trying to determine whether my shapefile that I pass is in M-aware (has M value). I keep getting an error

Attribute Error: 'Polyline' object has no attribute 'hasM'

Not sure what I am doing wrong.

import arcpy

def ValidInputGeometry(feat):
    is_valid_poly_line = False
    is_m_enabled = False
    error_msg = None
    if feat.type == "polyline":
        if is_m_enabled == feat.hasM:
            is_valid_poly_line = True
        else:
            error_msg = "The input geometry is not M-Aware"
    else:
        error_message = "The input geometry is not a polyline"
    return (is_valid_poly_line, error_msg)

station_interval = 10000
in_polyline_m_layer = "C:/Users/Algorithm/LWOB_levee.shp"
target_feature_rows = arcpy.SearchCursor(in_polyline_m_layer, "FID = 0", "", "", "")
layer_description = arcpy.Describe(in_polyline_m_layer)
shape_field = layer_description.ShapeFieldName

for row in target_feature_rows:
    feat = row.getValue(shape_field)

    (Valid, Msg) = ValidInputGeometry(feat)


texts: How to consecutively number polygons from north to south following a polyline?

texts:
I have a shapefile (polygon) with a large number of polygons (roughly 10,000) that form a narrow but very long belt that windes through half a country and another shapefile (polyline) that represents the western border of the belt.
Of those 10,000 polygons only about one third touches the border. So simply cutting the line into smaller lines using intersect, number those and then join back does not work.
I need to gibe each polygon a unique ID-number that consecutively increases from the first polygon in the north to the last polygon in the south.
My first idea, to simply create centroids, calculate the y-value for them and sort decreasingly, does not work either because the southern end is not the most southern point of the belt.

My question now is how to get this done with an automated action and without leaving the polygons behind that do not touch the border.

Edit:


texts:
We are using OpenLayers on mobile devices. When ever the keyboard is shown, OpenLayers resizes the map. This gives us performance issues when there are many features on that map. It takes a long time for the keyboard to appear/disappear.

Is there a way to prevent OpenLayers from resizing the map? Or are there any other ideas how avoid performance issues?

EDIT:
A user can click on a feature and an overlay is opened, where he can edit some of the values. The map is only in the background, mostly not visible for the user. So, the problem is not so much that the map gets bigger/smaller, but that the map is re-rendered / reloaded. So every time the keyboard is shown, the whole map is re-rendered (although invisible). Is there a way to stop re-rendering / reloading?


texts: Do outputs (shp, dwg) inherit vertical datum from input (las) files?

texts: I'm working in ArcMap. I pulled LiDAR data off earth explorer and created a DEM with the LiDAR data. I then created contours as a SHP file and exported to DWG (CAD). I was asked if the CAD file's vertical datum was NAVD88. I assume it is because I didn't change anything. If the LiDAR data was originally NAVD88, will the CAD file be in NAVD88?


texts:
The non-contiguous cartogram seems to solve a lot of the area distortions seen with traditional cartograms, while preserving the ability to identify areas by their shape lost with Dorling or Demers cartograms.

However, I can see how there might be problems related to the size of objects.  Are there any practical guidelines on when the ability to comprehend these cartograms breaks down, either by:

The ratio of the smallest to the largest value, or
The ratio of the smallest original area to the largest (e.g. in a map of US states, the viewer might have difficulty figuring out whether the values were encoded in the object size or rather in the object's size relative to its original size)

Or both?

Example cartogram with small states having large values:

Example cartogram with not much variation between values represented:


texts:
I have used QGIS2threejs for several years in QGIS2. Now we have upgraded to QGIS3 at work and the plugin no longer works.

The qgis2threejs window is blank and it displays the error “Any 3D objects not registered. There is a compatibility issue with QGIS 3D view. You need to close QGIS 3D view(s) and restart QGIS to use this preview.”

I have restarted QGIS and restarted my computer and tried to open QGIS2threejs as the first action. Still the same problem.

This page suggests going to view and manually closing any 3D map windows. I have checked this and I do not have any 3D map windows open.
https://github.com/minorua/Qgis2threejs/issues/147

I am using version 3.6.2. I am not able to change my version of QGIS because this is a work computer and I don’t have the permissions.

Any ideas how I can get the plugin to work?


texts: Pushing Geocoded Point to End of Parcel using ArcGIS Desktop?

texts:
I am trying to geocode for trees based on addresses. Basically I need the points to show up at the end of the parcel as opposed to the center (trees don't grow in the middle of houses). I have tried using the 'Side Offset' function in my ArcMap desktop (10.5) but this seems to confuse the program and a number of the points that were regularly geocoded disappear and nothing slides to the end of the parcel.

The geocoder I built uses building parcels and the 'General- Single Field' style. It works excellently until you try to offset anything. I have also tried this with street center lines (US Addresses- Dual Ranges) but this doesn't work as well. The image below shows in orange the parcel geocoded sites and the blue is the street geocoded sites.

The data I am trying to geocode has an address ex: "1501 N 15TH ST", and a zipcode ex: "19121" in two separate fields.

I am not opposed to writing code for this, but would like to see if there was an option to do this in the ArcMap geocoder tools before I went too far.


texts:
I have an array of Lines and I am using it to draw vectors in my map.
I want to replace two superposed Lines (or have short distance between them) with one Line.
Here is a picture that helps you to understand this problem:

The Input Lines :

The Output

PS: A Line is an ArrayList of points.


texts:
I am experimenting with CQL filtering in Geoserver. I have two layers, one (A) has 8000 features, the other one (B) has 3000 features or so.
I would like to search of layer A for intersection with layer B's features. Not all of Layer B's features, just a few based on their feature ID.

I've been searching around and I can't find a way yet to do the feature ID + geometry filter, I hope someone can suggest me how to do that.

So in the meantime, I have tried with another attribute I've put in, though, and the CQL Filter is currently like this with WMS:

layers=A
CQL_FILTER=INTERSECTS(the_geom,querySingle('B','the_geom','OBJECTID=6139'))

When I run this in my browser , the query goes forever and timed out. There is an error in the jetty : org.geoserver.platform.ServiceException: This requested used more time than allowed and has been forcefully stopped. Max rendering time is 60.0s

The example from Geoserver (the one using restricted layer vs bug sites) doesn't have this problem. I'm wondering if I'm writing the query wrong, whether my PC's just not up to it, or something else.  Should I be doing this in WFS instead?


texts: Adding and toggling legend with layer in OpenLayers?

texts:
I'm currently trying to add a legend to my map and to show/hide it, when the corresponding layer (an image layer) is shown/hidden.
My approach was to add a custom button with a .png of my legend. That worked without problems.

The image layer ("light") is deactivated when loading the map:

light.setVisibility(false);

Now, I want to show the legend when activating the layer and hide it when deactivating the layer and so on. And that works as well, but only once using these lines:

(function legend() {
    light.events.register("visibilitychanged", light, function() {
    if (light.getVisibility() === true) {
        map.addControl(button);
        } else {
        map.removeControl(button);
        }
    });
}());

I guess i have to reload the function somehow after the button is removed?

All things I tried didn't work and that was the best I was able to come up with.


texts:
I'm Working on a map and I want to put popups to show information about each queryable layer. This is what I've got so far, but I can't make it work.

HTML:

        X

Ol3:

/**
 * Elements that make up the popup.
 */
var container = document.getElementById('popup');
var content = document.getElementById('popup-content');
var closer = document.getElementById('popup-closer');

/**
 * Add a click handler to hide the popup.
 * @return {boolean} Don't follow the href.
 */
closer.onclick = function() {
  overlay.setPosition(undefined);
  closer.blur();
  return false;
};

/**
 * Create an overlay to anchor the popup to the map.
 */
var overlay = new ol.Overlay(/** @type {olx.OverlayOptions} */ ({
  element: container,
  autoPan: true,
  autoPanAnimation: {
    duration: 250
  }
}));

/**
* Layers of the map
*/

var wmsSource = new ol.source.ImageWMS({
  url: 'http://localhost:8080/geoserver/wms',
  params: {'LAYERS': 'ne:ne'},
  serverType: 'geoserver'
});

//Base layers

var IGAC = new ol.layer.Tile({
                type:'base',
                title:'IGAC',
                source: new ol.source.TileWMS({
                    url: 'http://geocarto.igac.gov.co/geoservicios/cien_mil/wms',
                    params: {LAYERS: 'Todas', FORMAT:'image/png'},
                    attributions:[
                        new ol.Attribution({
                            html: 'Instituto Geográfico Agustin Codazzi' +
                            'IGAC',
                        }),
                    ],

                })
            });

// Overlayers

var viviendas = new ol.layer.Tile({
                title:'Viviendas',
                source: new ol.source.TileWMS({
                    url:'http://localhost:8080/geoserver/FrCS/wms',
                    params:{LAYERS:'viviendas', format:"image/png"},
                })
            });

/**
* Create the map
*/

var vista= new ol.View({
                projection: "EPSG:3857",
                center: ol.proj.transform([-71.94012,4.3937],'EPSG:4326', 'EPSG:3857'),
                zoom:12.5,
            });

            var map= new ol.Map({
                target:'map',
                view: vista,
                layers:[ IGAC, viviendas],
                logo: {
                    href:"http://fundacionpervivir.com",
                    src:"img/pervivir_fondo.png",
                },
                overlays: [overlay],
                controls: ol.control.defaults().extend([
                    new ol.control.MousePosition({
                        coordinateFormat: ol.coordinate.toStringHDMS,
                    }),
                    new ol.control.ScaleLine (),
                    new ol.control.ZoomSlider(),
                ]),

            });

/**
 * Add a click handler to the map to render the popup.
 */
map.on('singleclick', function(evt) {
  var coordinate = evt.coordinate;
  document.getElementById('info').innerHTML = '';
  var vistaResolution = /** @type {number} */ (vista.getResolution());
  var url = wmsSource.getGetFeatureInfoUrl(
      evt.coordinate, vistaResolution, 'EPSG:3857',
      {'INFO_FORMAT': 'text/html'});

  if (url) {
    content.getElementById('info').innerHTML =
        '' + url + '';
  }
  overlay.setPosition(coordinate);
});


texts:
I am just trying to rename my file names from the list shp.import. But file.rename function doesn't helps.

shp.import <- list.files()
file.rename(from=shp.import, to=sub(pattern="dat_SPDF_", replacement="", shp.import))
file.rename(from=shp.import, to=sub(pattern=".shp", replacement="", shp.import))
shp.import

Here I need to remove 'dat_SPDF_' and '.shp' from the file name. I also didn't find any examples with stringr library.


texts: Sync a file geodatabase from a portable hard drive to mother database on my main computer

texts:
I have several file geodatabases to which I routinely add new data generated in the field.  To accommodate the mobile nature of this data creation and maintenance, I am working off of an external hard drive.  I would like the "mother" files of these geodatabases to live on my desktop, where I have a regular backup scheduled.  Is there a way I can set up a system to sync the updated geodatabases with this main GIS machine when I am back at the office?  I am unable to remotely access this main machine due to government restrictions.

I am running ArcGIS 10.2.2 with the ArcInfo license.

P.S. I am a new user and this is my first post, so I apologize if I'm missing anything or using poor form.


texts: Listing all exact Degrees, Minutes, Seconds for an entire area

texts: Does a database or dataset exist that lists all exact Degrees, Minutes, and Seconds (rounded to nearest decimal) for the continent of South America? What I am seeking is a list of DMS's that covers that entire geography; not just a few specific coordinates, but all of them. Given that each DMS represents ~60 nautical miles apart, the data shouldn't be overly complicated to form. The trouble is, I can't locate the consolidated data anywhere on the web!


texts: I am trying to access the fields which I can see in this particular layer. I am able to see that it has fields but fetching the values in those fields is something that I am unable to do. Is there someway to do that?


texts: Mapcache: ValueError: unsupported format character 'j' (0x6a) at index 55

texts:
For my first attempt at using Mapcache, I wish to create a caching proxy of Bing Maps, serving a WMS. I started 'mapproxy-util serve-develop' with the following .yaml configuration, assembled from studying the documentation :

services:
  wms:
    md:
      abstract: ''
      title: Microsoft Bing orbital imagery caching proxy
    srs: ['EPSG:3857']

layers:
-   name: Bing_imagery_layer
    title: WMS layer from Microsoft Bing orbital imagery tiles
    sources: [Bing_imagery_cache]

caches:
  Bing_imagery_cache:
    grids: [GLOBAL_WEBMERCATOR]
    sources: [Microsoft_Bing_imagery]
    cache:
      type: file
      directory: /home/jm/projects/mapproxy/cache
      directory_layout: tms

sources:
  Microsoft_Bing_imagery:
    type: tile
    grid: GLOBAL_WEBMERCATOR
    url: http://ecn.t2.tiles.virtualearth.net/tiles/a%(quadkey).jpeg?g=1&mkt=en-US

And here is the observed behaviour:

$ mapproxy-util serve-develop Bing_orbital_imagery_cache_WMS.yaml
[2016-05-24 00:17:50,048] mapproxy.config - INFO - reading: /home/jm/projects/mapproxy/mymapproxy/Bing_orbital_imagery_cache_WMS.yaml
[info]  * Running on http://127.0.0.1:8080/ (Press CTRL+C to quit)
[info]  * Restarting with reloader
[2016-05-24 00:17:50,771] mapproxy.config - INFO - reading: /home/jm/projects/mapproxy/mymapproxy/Bing_orbital_imagery_cache_WMS.yaml
[2016-05-24 00:17:54,457] mapproxy.wsgiapp - CRITICAL - fatal error in wms for /wms?FORMAT=image/jpeg&VERSION=1.1.1&SERVICE=WMS&REQUEST=GetMap&LAYERS=Bing_imagery_layer&STYLES=&SRS=EPSG:3857&WIDTH=512&HEIGHT=512&BBOX=-1232777.2003892,1692622.3427892,-1227885.2307892,1697514.3123892
Traceback (most recent call last):
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/wsgiapp.py", line 189, in __call__
    resp = self.handlers[handler_name].handle(req)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/ows.py", line 38, in handle
    return self.services[service].handle(req)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/base.py", line 30, in handle
    return handler(parsed_req)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/wms.py", line 131, in map
    renderer.render(merger)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/wms.py", line 530, in render
    layer_merger)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/wms.py", line 567, in _render_capture_source_errors
    raise ex[1]
ValueError: unsupported format character 'j' (0x6a) at index 55
Traceback (most recent call last):
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/wsgiapp.py", line 189, in __call__
    resp = self.handlers[handler_name].handle(req)
  File "/home/jm/projects/mapproxy/local/lib/python2.7/site-packages/mapproxy/service/ows.py", line 38, in handle

In the URL passed to Mapcache, the only occurrence of the letter "j" mentioned in the error message is in "FORMAT=image/jpeg" - so I thought maybe the forward slash must be encoded as %2F but trying that yielded the same behavior.

Even stranger, I tried passing "FORMAT=image/png" so that the entire URL contains no occurrence of the letter "j" and I got exactly the same result: "ValueError: unsupported format character 'j' (0x6a) at index 55"

At this point, I am at a loss for new hypothesis... From the lack of Google results for "unsupported format character" "mapcache" it seems that I managed to err in exotic ways..

If any of you has hypothesis about where I should look, I'll be thankful.


texts: Create lines between two points, given an ID

texts:
I have a spatialite point table (p1) with an additonal column (pid) that relates each point to a point in another point table (p0) in the same database. I want to create a line between every pair of points where p1.pid matches p0.id. This is a 1:1 relation.

How can I do that in spatialite or in QGIS 2.14?


texts: I have a KML file with a LineString. I'm trying to load the points in Python. I'd like to use the FastKml library, but am open to other solutions. FastKml has one example for reading KML, but I have not been able to convert that to get the LineString points.


texts: how to find angle between two points in QGIS?

texts:
I want to find angle between two points(LatLong) in QGIS. For example

Point A = 80.2065,13.0876
Point B = 80.2064,13.0874

My CRS is WGS84 EPSG : 4326.

I am trying to provide driving direction this requires angle between points.


texts: how to differentiate between common parts of the triangle and extract linestring from it

texts:
Can you help me break down shape A to different linestrings: lines 2 & 3 AND line 1?  The shape contains three different polygons (A, B, C). I wish to determine the start/end point of each part with the common points of the other polygons.

I'm using Postgis 1.5 / Python Shapely


texts: Which method is most suitable to change the pixeltype of the raster in PostGIS?

texts:
I have saved a raster in the PostGIS database and want to select it as PNG. However, ST_AsPNG only accepts rasters with pixel type 8 BUI or 16 BUI. But in this case it is a raster with pixel type 32 BF.
The only possibility I have seen so far is ST_Reclass. But you need to define a 'reclassarg' which makes it unnecessarily complicating for my usage, because I only want to change the pixel type.
Does anyone know a more suitable method for this purpose?


texts:
am trying to create a map, which shall show NRW ( a part of Germany) with some extra gimmiks. Those Gimmiks are Coordinate based locations of several Manufacturers. Later on i want to adjust the background color to statistik Data.

The Base map is Ready.

Now i tried to add the Locations of the manufacturers. Therefore i created a Textdata with the coordinates. via layers/Add Layers / Import Textdata as Layer (I have a german version so it is free translated)i tried to add this Layer. After several tries it worked
The Problem is now, that this Layer is not overlapping with the Card i have so far created. I think it has something to do with the reference Coordinate System but i dont know how to correct it.

Can somebody pls help me?


texts: How to remove features from a Leaflet GeoJSON layer?

texts:
In a Leaflet map, it is possible to add features to a GeoJSON layer using the addData method:

  L.geoJSON: Creates a GeoJSON layer. Optionally accepts an object in
  GeoJSON format (you can alternatively add it later with addData
  method)

  addData(data): Adds a GeoJSON object to the layer.

Is it possible to remove all features from the GeoJSON layer, without destroying the layer and recreating it?


texts: How to use "Data-Driven Input Mask" QGIS plugin?

texts: I have installed Data-Driven Input Mask plugin in QGIS and added necessary relationships between PostgreSql tables. But I don't find any changes in the UI. How to enable/configure this plugin functionality in the attribute edit form as mentioned in the documentation?


texts: Creating a 5x5 poylgon in ArcGIS with a predetermined center with coordinates X, Y

texts:
I need to build a polygon (5x5 km) centered on a coordinate X,Y. I want to use this 5x5 square for an eventual clip.

How should I proceed in ArcMap?


texts: how to calculate area if it is falling in two utm zones

texts: I am using ArcGIS 10.1. How do I calculate an accurate area measure for a polygon, where portions of it are falling in two zones, i.e. UTM ZONE 44 & 45 (India)?


texts: How to convert my GML/KML/SHP files into deegree workspace?

texts:
I am newbie in geographic and want to convert my GML file to deegree workspace.

What I did so far

I have download deegree and install it on local.
Run the application and its service

Now, I want to convert my GML/KML/SHP files to Deegree workspace.

Please help me or any suggestion what i do.


texts: Clip multiple line layers with several polygon layers at the same time

texts:
I am in the process of making the database for 9 consecutive years in 150 local administrative units. To do that I have to clip 9 line layers (roads for 9 years) with 150 polygon layers.

So far I have been doing it with the clipping tool in Batch process, but there are several problem I would like to solve a bit more easier.

Problem 1) When I clip the line layers I have to manually update the length of the lines. Since there are 1350 cases (9 years x 150 administrative units), it would take me a 2-3 days of repetitive and boring work that I would like to avoid. (Question 1) So the question is if there is any tool/plugin that can clip multiple line layers with multiple polygon layers and create separate new line layers for every year and administrative area. (Question 2) Or, is there any tool/plugin to update all the new layers at the same time

Problem 2) When I am done with clipping and updating of the line lengths, (Question 3) is there any tool/plugin that allows me to export all the line layers into one excel or CSV file?

I am using QGIS Desktop 3.10.2 with GRASS 7.8.2


texts: Save Google Map layer as in QGIS using Python script

texts:
I am trying to save Google Maps by using save raster layer as via Python script:

urlWithParams = 'type=xyz&url=https://mt1.google.com/vt/lyrs%3Ds%26x%3D%7Bx%7D%26y%3D%7By%7D%26z%3D%7Bz%7D&zmax=19&zmin=0'
rlayer = QgsRasterLayer(urlWithParams, 'Google Satellite', 'wms')
renderer = rlayer.renderer()
provider = rlayer.dataProvider()
pipe = QgsRasterPipe()
pipe.set(provider.clone())
pipe.set(renderer.clone())
extent = self.iface.mapCanvas().extent()
crs = self.dockwidget.crs_select.crs()
file_writer = QgsRasterFileWriter(self.dockwidget.input_path.text())
file_writer.Mode(1)
file_writer.writeRaster(pipe, provider.xSize(), provider.ySize(), extent, crs)

It works, but the result not as correctly as I do it by Export in QGIS.

Any ideas?


texts: Getting "Workspace Already in Transaction Mode" Error when trying to use updateCursor on versioned sde

texts:
There are are numerous post on this error, but I can't seem to find a solution for myself. We have a versioned sde and need to add new fields to some of our feature classes and populate them with an updateCursor.

code here:

import arcpy
print('Script Started')

arcpy.env.workspace = r"Database Connections\Connection to gisprd.sde"
edit = arcpy.da.Editor(arcpy.env.workspace)
edit.startEditing(True, False)

fcNames = [
    "Transformer",
    "DynamicProtectiveDevice",
    "Switch",
    "BulkDeliveryPoint",
    "BusBar",
    "CableFault",
    "ConduitSystem",
    "Fuse",
    "NetworkProtector",
    "OpenPoint",
    "OverheadPrimaryLine",
    "OverheadSecondaryLine",
    "Riser",
    "ServicePoint",
    "SupportStructure",
    "SurfaceStructure",
    "UndergroundPrimaryLine",
    "UndergroundSecondaryLine",
    "UndergroundStructure",
    "VoltageRegulator",
    "WireChange"
    ]

fcs = arcpy.ListFeatureClasses("","","LIPA.Distribution")
print('Starting loop')
for fc in fcs:
    # get feature class name
    name = arcpy.Describe(fc).baseName[5:]
    print(name)
    path = arcpy.Describe(fc).catalogPath
    # only do to the features listed above
    if name in fcNames:
        # Add field to store static object ids
        arcpy.AddField_management(path,'ORIG_OID', 'LONG')
        oid = [] # initalize empty to store objectids
        idx = 0 # initalzie empty index counter to add correct objectid to new field for static oid
        # Use cursors to get actual objectids for features of interest
        print('Starting Search')
        with arcpy.da.SearchCursor(path, ('OBJECTID',)) as cursor:
            for row in cursor:
                oid.append(row[0])
        print('Search finished, starting update')
        # Iterate update cursor over new field
        edit.startOperation()
        with arcpy.da.UpdateCursor(path, ('ORIG_OID',)) as upCursor:
            # for each row update with the correct object id and then increment the counter index
            for row in upCursor:
                row[0] = oid[idx]
                idx += 1
                upCursor.updateRow(row)
        print('Finished update')
        edit.stopOperation()
        del cursor
        del upCursor

edit.stopEditing(True)
del edit

I got the following traceback:

for row in upCursor:
RuntimeError: workspace already in transaction mode

Can anyone assist?

links of similar questions:

Debugging RuntimeError: workspace already in transaction mode from arcpy.da.UpdateCursor and ArcSDE feature classes?

Workspace already in transaction mode

The above two are both trying to work in different versions. We only do stuff like this on weekends after all the work is done for the week, everything has been reconciled and compressed (database is at state 0; not worried about versions).

Checking via ArcPy if ArcMap is in edit session?

This one seems to indicate that error means an edit session is not open? So maybe I am setting up the editor wrong?


texts: From hundreds of shapefiles to PostGIS with some kind of directory/file structure

texts:
We have more than 700 shapefiles which we will import in a PostgreSQL/PostGIS database.
We want to use a single database to make joins on our tables without using foreign database wrapper.
Our data guy uses directories to classify shapefiles for different theme (water, energy...) so he is not comfortable with the idea of having all data in one database.

In PostgreSQL the only way I see is to use schemas but we'll have only one level of depth.

How do you handle so many tables in one database without losing too much time looking for it in your everyday routine? Do you know a tool to manage a catalog of tables in a "directory/file structure" way? How would you change your workflow from a lot of shapefiles to one big database?


texts: refine nearest_node with radius-distance ( ST_Distance, ST_DWithin )

texts:
i am a hobby-IT-person and currently stuck with some simple problem, maybe somebody can help me:

This is my working POSTGIS-function, source by underdark:

ALTER TABLE distance_graz
   ADD COLUMN nearest_node integer;

CREATE TABLE temp AS
 SELECT a.osm_id, b.id, min(a.dist)
 FROM
(SELECT planet_osm_point.osm_id, min(ST_Distance(planet_osm_point.way, distance_graz.the_geom))  AS dist
      FROM planet_osm_point, distance_graz
      GROUP BY planet_osm_point.osm_id) AS a,

     (SELECT planet_osm_point.osm_id, distance_graz.id,
             ST_Distance(planet_osm_point.way, distance_graz.the_geom) AS dist
      FROM planet_osm_point, distance_graz) AS b
   WHERE a.dist = b. dist
         AND a.osm_id = b.id
   GROUP BY a.osm_id, b.id;

UPDATE distance_graz
SET nearest_node =
   (SELECT id
    FROM temp
    WHERE temp.gid = distance_graz.gid);

What it does is to get the node "a" with the minimal distance to another node "b" and write it in a new table

But now I want to refine this by only getting "nearest points "b"" if they are within a certain distance (5000m) to this point "a".

I think the best way to do this is by adding "ST_dwithin" to the quer:

SELECT planet_osm_point.osm_id, ST_DWithin(planet_osm_point.way,distance_graz.the_geom,5000) AS wit

The Question now is how and where do i implement this?


texts: Did Digitizing tools last update remove polygon cutting?

texts:
In a previous project i used the Digitizing Tools plugin to cut polygons from other polygons.
But it seems that with a new update that feature is removed? I can't seem to use the most of the buttons.
Any advise on a different plugin that will easily help me cut polygons from other polygons?
Especially when its difficult with underlying polygons that get completely covered by the polygon on top, in the same layer.
I am not a professional with QGis so I hope there is an easy fix for this.

To specify; I use QGIS Desktop 3.8.0 (Zanzibar) with GRASS 7.6.1.
The Digitizing Tools plugin is updated to the latest version 1.5.1

Here is an image to shhow how i can't use the tool dispite being in editing mode for the layer.


texts: How do I get a decent easy to use watermark from an image in Quantum?

texts: Our users need to be able to put a graphic watermark onto maps they create in the print composer. Currently we get them to georeference the image as a raster and then create a transparency. This is quite complex for some of my beginners. Anybody got any suggestions?


texts: How to get different clazz/flags for roundabout?

texts:
I'm using osm2po to import Open Street Map data. I want to get information about roundabouts. In the final SQL file I need to know if a road is part of a roundabout or not.

The only way to get this information in the SQL file is to change the clazz or flag value.

How can I get a different clazz value for roundabouts?

I have tried to change the configuration file but without any effect :

wtr.tag.junction.roundabout =    1, 77, 35

I never get 77 value in clazz column.


texts:
I'm starting to work with the advanced findlabel function in ArcMap label expression editor. I will say I'm a complete novice at this up front. Here is what I have so far:

Function FindLabel ( [SwineWeFin] , [SwineGrFin] , [SwineNurse] , [SwineGesBo] , [SwineSowLt] , [SwineGilti] , [CattleDaHl] , [CattleDaJr] , [CattleDaHf] , [CattleBeFn] , [CattleBeFe] , [CattleVeal] , [TurkeyFin] , [TurkeyPlBr] , [ChickenLBr] , [ChickenPul]     )

 if ( [SwineWeFin]  > 0) then  FindLabel = [SwineWeFin] & " WTM "

 if ( [SwineGrFin] > 0) then  FindLabel = [SwineGrFin] & " FTM"

 if ( [SwineNurse]  > 0) then  FindLabel = [SwineNurse] & "Nursery"

 if ( [SwineGesBo] + [SwineSowLt]  > 0) then  FindLabel = [SwineGesBo] + [SwineSowLt]  & "Sows"

 if ( [SwineGilti]  > 0) then  FindLabel = [SwineGilti]  & " Gilts"

 if ( [CattleDaHl] +  [CattleDaJr] +  [CattleDaHf] +  [CattleBeFn] +  [CattleBeFe] +  [CattleVeal]  > 0) then FindLabel = "Cows"

 if ( [TurkeyFin] + [TurkeyPlBr]  > 0) then  FindLabel =  "Turkey"

 if ( [ChickenLBr] + [ChickenPul]  > 0) then  FindLabel = "Chickens"

End Function

This mess actually works for me. The problem that I have is when two or more of the "if' criteria apply to the same point.  Currently, the label will only show the first criteria that is found to be true.  I want all criteria that are true to be shown as a label.  I would prefer if it showed each TRUE item on a new line for the label.

Is what I'm wanting to do even possible? And if so, how would I script it? Thanks!!!!


texts: Determine on what side of a line my current position is at, no matter orientation

texts:
I know that this question has been up before, but I'll give it a shot since it's not exactly the same question.

Imagine a straight line consisting of a start point and an end point.
The orientation of the line in 2D can be whatever we choose.
Now imagine that we are traveling along this line, but not exactly on the line.
I would like to give the directional (e.g right, left) feedback so that they know what direction to head to stay on the line.

I've been reading and trying this solution:
https://math.stackexchange.com/questions/274712/calculate-on-which-side-of-straign-line-is-dot-located

I'm currently stuck when trying to grasp what happens if the we travel in the opposite direction.. I's there a solution that will always give us the correct direction?
No matter the line's orientation and our traveling direction.
I have the possibility to place one or more reference points in the plane that are known.


texts:
I'm trying to orthorectify an image with the following commands:
gdalwarp -rpc -to RPC_DEM=dem.tif input.tif output.tif
gdalwarp -rpc input.tif ortho_output.tif
I'd expect the "ortho_output.tif" to look different, especially in non-flat areas, but both images are identical. The RPCs in the metadata seem to work because both images are warped to the correct locations, but the RPC_DEM option doesn't seem to do anything.
Am I missing something? The processing time is slightly longer when using the DEM, so I assume something is happening. I've also experimented with using RPC_HEIGHT=500, but to no affect.
My gdal version is 3.0.4, I've also checked that my DEM (in this case SRTM) covers the full extent of my image and that it is using ellipsoidal heights.
Are there any public datasets I could use as test data?


texts: Select features from another layer based on a selection in PyQGIS

texts:
I have two polygon features: epoca and result.

I would like to select the features that correspond to the number 1 feature (features A and B) and then select the corresponding features with the number 2 feature (features C and D).

I had done the following:

from processing.core.Processing import Processing
import processing.tools

layerResult = processing.getObject('lyr_result')
layerEpoca = processing.getObject('lyr_epoca')

processing.tools.general.runalg('qgis:selectbylocation', layerResult, layerEpoca, u'equals', 0.005,0)

layerResult.invertSelection() #I'm searching for the different polygons
selectFeatures = layerResult.selectedFeatures() #are the features selected in result

processing.tools.general.runalg('qgis:selectbylocation',layerEpoca, layerResult, u'within', 0.005,0) #here select the features in epoca all at once!

What is being done is that you are selecting features A, B and C, D at the same time. I would like to select these features at different times. Is this possible?

I had thought of something like this:

for feat in selectFeatures:
    processing.tools.general.runalg('qgis:selectbylocation',layerEpoca, selectFeatures, u'within', 0.005,0)
#----but it did not work----#


texts: Is there machine learning for image classification in ArcGIS 10.3?

texts:
I have access to ArcGIS 10.3. Do I need to get ArcGIS Pro to have access to machine learning?
I've read that this is a built-in function in ArcGIS Pro.
links:
http://www.esri.com/esri-news/arcuser/spring-2018/machine-learning-in-arcgis
https://www.youtube.com/watch?v=duZ7jor_YrU


texts: Print Composer Template produces empty map?

texts:
I have a print composer template which I'm trying to populate with different subset queries from a CSV-derived vector layer using the Python API. However, all the output PDFs contain everything except the layers on the map! I'm doing this with a standalone app, and it functions fine in that it does produce PDFs and the data is all loaded and interrogated just fine (I can subset query, count features etc.)

Here's the method I'm using to render the composition, as you can see I've tried to set the layers on the map in different places (actually there are 2 maps - one shows the extent rectangle of the other, but I can't get either to show anything) - this is my first QGIS scripting experience so I'm sure I'm missing something very obvious...

def compose(layers, species_binomial):
    # vlayer is the CSV vector layer, outline is the ADM0 country vector.
    vlayer, outline = layers

    # construct a QStringList for the renderer setLayerSet
    layers = QStringList(vlayer.id())
    layers.append(outline.id())

    # 1 maybe set the registry layers???
    QgsMapLayerRegistry.instance().addMapLayer(vlayer)
    QgsMapLayerRegistry.instance().addMapLayer(outline)

    # 2 maybe set the canvas layer???
    canvas = QgsMapCanvas()
    canvas.setExtent(outline.extent())
    canvas.setLayerSet([QgsMapCanvasLayer(outline), QgsMapCanvasLayer(vlayer)])

    # 3 set the renderer layers???
    renderer = canvas.mapRenderer()
    renderer.setLayerSet(layers)
    renderer.setExtent(outline.extent())

    composition = QgsComposition(renderer)
    with open(os.path.join(ROOT, 'species_composer.qpt'), 'rt') as fp:
        composer_template = fp.read()
    doc = QDomDocument()
    doc.setContent(composer_template)
    composition.loadFromTemplate(doc)

    # 4 maybe set the map layers of one of the maps???
    map1 = composition.composerMapItems()[1]
    map1.setLayerSet(layers)

    title = composition.getComposerItemById('title')
    title.setText(" ".join(species_binomial))  # this works fiiiine

    out_path = os.path.join(ROOT, 'out',  "%s_%s.pdf" % species_binomial)

    # make a lovely blank PDF with a title and a blue map rect
    ...

What am I doing wrong? Here's an example output...


texts: Using QGIS Zonal Stats Plugin from Python Console?

texts: I am using the Zonal Stats plugin in QGIS to extract raster statistics from overlaying polygons (I have shapefiles of species distribution and I want to extract environmental data from within each species' range). I have 300 or so files I need to get data from and so would like to write a script to run in the python console, however I am a complete novice with python and have no idea how to do this.


texts: FeatureClass error message when running ArcGIS Python tool

texts:
I am trying to remove rows from a feature class that is included in a python tool in ArcGIS 10.7. I have run the code as an independent Python script and it runs fine, but when I run the following code, I get an error. tmpfile is a copy of an existing feature that is in a geodatabase.
tmpfile = os.path.basename(tempfile.TemporaryFile().name)
arcpy.CopyFeatures_management(infile, tmpfile)

### remove unnecessary values
vals = ["Unknown",""]
for val in vals:
    query = """ "DATEDISC" = '%s' """ %val
    tmpfile = arcpy.SelectLayerByAttribute_management(tmpfile,"NEW_SELECTION", query)
    arcpy.DeleteRows_management(tmpfile)

In the standalone script, this removes any rows with the values of "Unknown" or "". In the tool, I receive the error message for arcpy.SelectLayerByAttribute_management() that
The value cannot be a feature class
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

This seems like odd behavior for this function.
Does anyone know why this would be happening?


texts: Creating random points within polygon from database?

texts:
Although I have a solution for it in ArcGIS named Hawth's tools, but I wanted to know whether it is possible in QGIS or not. Additionally, I am not having a lower version of ArcGIS on which I can do my work:

I have a database file having village names
I have village shapefile having names
I want to geocode the database file on the basis of matching names for villages, but I want to geocode it randomly in the polygon and not overlapped on each other where names are matching.


texts: What does it mean when I get rows executed, but with no resulted table?

texts:
While applying the following tutorial: http://anitagraser.com/2011/02/09/creating-catchment-areas-with-pgrouting-and-qgis/ utilizing vector network which has been tested (and demonstrated functionality) for shortest path application, it returned the following message:

  Query returned successfully: 1013 rows affected, 199 ms execution time

but returned no table as output!

The code I used is:

create table catchment_1 as
select
    id,
    the_geom,
    (select sum(cost) from (
       SELECT * FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(geom) as cost FROM schema.cvl',
    1, 320, false, false
  )) as foo ) as cost
from schema.cvl_vertices_pgr

I also tried passing "id" to the function instead of "320" as per the tutorial but got the following msg:

  ERROR: function pgr_dijkstra(unknown, integer, bigint, boolean, boolean) does not exist SQL state: 42883

Although, pgr_dijkstra is working well with other codes.. so, it's one of the recognized functions in my installed pgRouting layer

I'm using QGIS 2.0.1 and pgRouting pg93-binaries-2.4.0


texts: how to create temperature heat map on a web?

texts:
I'm work on a website of weather data .

I need to make heat map with temperature data to show the status of temperature in the state

my data(lat,lng,temp), I'm trying to work with Google maps but it embed only  (lat,lng) .

can any one know how to make it or know other applications can help me ?


texts: Connecting CartoDB to existing MySQL database?

texts: I'm trying to create a simple heat map based on longitudinal and latitudinal map coordinates that can be transferred from a database on MySQL that I've created. I read that CartoDB can be connected to other databases, but I'm not sure how to go about doing that.


texts: Rendering using color ramp for negative value in PyQGIS

texts:
I'll try to rendering single band raster (GEBCO bathymetry) by using PyQGIS..
This is my script so far

bathymetry_file = "E:/byth_gebco.tif"

layer_bathymetry = QgsRasterLayer(bathymetry_file, 'Raster Byth')

s = QgsRasterShader()
c = QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
c.setClassificationMode(CONTINUOS)
c.setMaximumValue(625)
c.setMinimumValue(-6244)
i = [
    QgsColorRampShader.ColorRampItem(0, QtGui.QColor('#f7fbff'), '0'),
    QgsColorRampShader.ColorRampItem(-200, QtGui.QColor('#deebf7'), '-200'),
    QgsColorRampShader.ColorRampItem(-1000, QtGui.QColor('#c6dbef'), '-1000'),
    QgsColorRampShader.ColorRampItem(-2000, QtGui.QColor('#9ecae1'), '-2000'),
    QgsColorRampShader.ColorRampItem(-3000, QtGui.QColor('#6baed6'), '-3000'),
    QgsColorRampShader.ColorRampItem(-4000, QtGui.QColor('#4292c6'), '-4000'),
    QgsColorRampShader.ColorRampItem(-5000, QtGui.QColor('#2171b5'), '-5000'),
    QgsColorRampShader.ColorRampItem(-6000, QtGui.QColor('#08306b'), '-6000')
]
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps = QgsSingleBandPseudoColorRenderer(layer_bathymetry.dataProvider(), 1, s)
layer_bathymetry.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(layer_bathymetry)

This is the expected result:

But this is what i've got from that script

It seems like QGIS API can't handle the negative pixel value,
what should i do to render that bathymetry remind that bathymetry has negative elevation?


texts: Joining attributes from two different polygon layers

texts:
I'm trying to join the attributes from the polygon layers shown on the images.

The result I'm looking for is one layer that contains the square and the attributes of both triangles: the triangles are perfectly overlapped with the square.

I've tried to use "Join attributes by location" but the output layer is not what I'm looking for as it contains two overlapped squares one for each triangle.

Also, I've tried to intersect the layers but the result are two triangles and to join the vectorial layers but the result is again two triangles.


texts: Error after adding extensions to GeoServer

texts:
I had installed java8, tomcat8.5.12, geoserver 2.9.4 in CentOS 6.9. I had added gdal libraries to the GeoServer from http://docs.geoserver.org/stable/en/user/data/raster/gdal.html.

I had restarted the tomcat there is no additional extensions added to the stores. Also, when I checked logs, I could see this:

it.geosolutions.imageio.gdalframework.gdal utilities load gdal failed to load the GDAL native libs.
this is not a problem unless you need to use the GDAL plugins:they wont be enabled.

I think it seems to be some path configuration error.


texts: Getting PNG binary data from GDALDataset

texts: I have an in memory GDALDataset of driver type PNG. How can I get the binary PNG data to send to a web client without first writing the GDALDataset PNG to disk and reading it back as a binary file?


texts: Which spatial weighting function does ArcGIS' GWR use?

texts:
I am running Geographically Weighted Regression and originally was using ArcGIS to carry out the analysis until I read/heard mixed reviews about it. I have now tried a different programme to run it (SAM) and get different results.

Does anyone know what the spatial weighting function is in ArcGIS? I think it is Gaussian but am not 100% sure. I cannot select this with Arc but there are different options to chose from with SAM.

Thanks!


texts:
I'm developing a Web Application that is using the ArcGIS Javascript API.  When adding a Feature Layer I'm using the setDefinitionExpression to limit the features that are returned.  The expression used can be quite large (greater than 2000 in length).
Because the expression can be large, I am attempting to use the Proxy page, however I am seeing unusual results.  I am finding that the proxy page is not always used, even though I am using the exact same query expression.  I have conformed this by utilizing IE's debugger tool and also using Fiddler.

In my application if I add 2 layers (Tiled and Dynamic) prior to adding the Feature Layer I consistently see that the Proxy page is not being used. I can see via the debugging tools that a POST is being performed directly to my REST service.  The content-length is over 6000.

If I add just the Tiled layer prior to adding the Feature Layer, I can see that an attempt is made to use the Proxy page.  The proxy page is having issues, therefore my Feature Layer does not get added.

Can someone please explain why the Proxy page is not always used even though it would seem that it should be, given the length of the query string used in the setDefinitionExpession.  I am very confused as to why layers added prior to adding a Feature Layer would affect the use of the proxy and why adding the Feature Layer doesn't fail when the proxy is not used.

The first function below shows where I am setting the esriConfig.defaults.
The second section of code is from my class (tpsMap) where I am adding a Feature layer.

require(["tps/map",
    "esri/config",
    "dojo/domReady!"],
    function (tpsMap, esriConfig) {

        // Setup Esri Config
        esriConfig.defaults.io.proxyUrl = "http://phunt-pc/EsriProxy/proxy.ashx"
        esriConfig.defaults.io.alwaysUseProxy = false;

        // Create and display Map
        tpsMap = new tpsMap(oMapConfiguration);
        tpsMap.mapQuery = oMapQuery;
        tpsMap.createMap("divMap");
    });

this._featureLayer = new FeatureLayer(this.mapQuery.FeatureLayer.url,
       {

            mode: FeatureLayer.MODE_SNAPSHOT,
            opacity: 0.50,
            infoTemplate:    this._getInfoTemplate(this.mapQuery.FeatureLayer.InfoTemplate),
            outFields: this.mapQuery.FeatureLayer.OutputFields,
            className: "feature"
        });

       // Limit the features based on the expession
  this._featureLayer.setDefinitionExpression(this._getQueryExpression());

  this.map.addLayer(this._featureLayer);


texts:
I have a .tif file and .ers .tfw with the same name, I would like to transform it to a GeoTIFF to raster it using gdal2tiles.
the .ers file content looks like:
DatasetHeader Begin
    Version     = "7.0"
    Name        = "705.ers"
    SourceDataset   = "705.tif"
    LastUpdated = Tue Jul 17 18:12:26 GMT 2012
    DataFile    = "705.tif"
    DataSetType = Translated
    DataType    = Raster
    ByteOrder   = MSBFirst
    CoordinateSpace Begin
        Datum       = "SIRGAS"
        Projection  = "SUTM24"
        CoordinateType  = EN
        Rotation    = 0:0:0.0
    CoordinateSpace End
    RasterInfo Begin
        CellType    = Unsigned8BitInteger
        CellInfo Begin
            Xdimension  = 0.080000000000000002
            Ydimension  = 0.080000000000000002
        CellInfo End
        NrOfLines   = 7425
        NrOfCellsPerLine    = 7638
        RegistrationCoord Begin
            Eastings    = 361854.52000000002
            Northings   = 7829895.6399999997
        RegistrationCoord End
        RegistrationCellX   = 0.5
        RegistrationCellY   = 0.5
        NrOfBands   = 3
        BandId Begin
            Value       = "Red"
        BandId End
        BandId Begin
            Value       = "Green"
        BandId End
        BandId Begin
            Value       = "Blue"
        BandId End
        RegionInfo Begin
            Type        = Polygon
            RegionName  = "All"
            SubRegion   = {
                0   0
                0   7425
                7638    7425
                7638    0
            }
        RegionInfo End
    RasterInfo End
DatasetHeader End

and .tfw file content are:
0.080000000000000002
-0
-0
-0.080000000000000002
361854.52000000002
7829895.6399999997

Well, the .tif is 170MB so, not possible to upload or show it here.


texts: ArcPy - Get Date as parameter value and calculate Date field

texts:
In my ArcPy script I have an input parameter for the Date, so the user is asked to input a date value (time is not needed) using the data type GPDate.
In the script I calculate a new field, field data type DATE, which should include for all features the date that the user entered.
I cannot find a way how to populate this new field with the value of the input parameter.
My code:
date = parameters[0].valueAsText
arcpy.management.CalculateField(myFC, "myField", date, "PYTHON3", "", "DATE")

I already tried 100 varieties with datetime, strftime and different formatting of the date (my default input is dd.mm.YYYY), but either I get an error indicating "unexpected EOF while parsing" or that the value cannot be written to the attribute table.


texts: GIS Application for online map and journey planner

texts:
I am planning to implement an online map application with journey planner for a city. I don't have any experience in GIS development but I have a few years of experience in Java development.

So could you guys suggest me how I could start off my implementation. Any recommendation for books or documentations? What kind of technology stack should I use?


texts: plotting a SpatialPolygons object using lines() without lines connecting separate polygons

texts:
This is almost an extension of a question asked on stackoverflow (https://stackoverflow.com/questions/12196440/extract-feature-coordinates-from-spatialpolygons-and-other-sp-classes), so I'll use the demo dataset from there.

I'm attempting to plot these polygons with lines(), so I can overlay these with an image.plot. To do this, I use fortify, like so

    # taken from the linked stackoverflow question
    Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
    Sr2 = Polygon(cbind(c(5,4,2,5),c(2,3,2,2)))
    Sr3 = Polygon(cbind(c(4,4,5,10,4),c(5,3,2,5,5)))
    Sr4 = Polygon(cbind(c(5,6,6,5,5),c(4,4,3,3,4)), hole = TRUE)

    Srs1 = Polygons(list(Sr1), "s1")
    Srs2 = Polygons(list(Sr2), "s2")
    Srs3 = Polygons(list(Sr3, Sr4), "s3/4")
    SpP = SpatialPolygons(list(Srs1,Srs2,Srs3), 1:3)

    image.plot(x = c(0:10), y = c(0:10), z = matrix(runif(100, 0,1), nrow = 10),
        col = terrain.colors(20)) # fake data so lines() plays nice.
    lines(fortify(SpP))
    # there's a weird line connecting the square island in the middle to the polygon surrounding it

What I would like to get is an image like this, which I can produce with ggplot...

    require(ggplot2)
    ggplot(aes(x = long, y = lat, group = group), data = fortify(SpP)) + geom_path()

I can see two possible solutions to this, either an analog to the "group =" argument in the aes() call in ggplot(), or some way to get fortify to insert 'NA' between different groups (this is what the lines() function interprets as breaks). I can produce an effect manually, but I couldn't do this with the real dataset.

    SpP3 <- rbind(SpP2[1:5,], rep(NA, 7), SpP2[6:10,], rep(NA, 7), SpP2[10:14,], rep(NA, 7), SpP2[15:19,])
    image.plot(x = c(0:10), y = c(0:10), z = matrix(runif(100, 0,1), nrow = 10), col = terrain.colors(20))
    lines(SpP3)

edit: I'd like to be able to produce these plots using image.plot, since it's faster for my real dataset, and I will be more compatible with what others I'm working with are using.

Thanks!


texts: Iterating over table to create chainages using Python in ModelBuilder

texts:
I routinely have to do the following workflow:

Map GPS points;-
Create lines based on those GPS points;
Split the lines at a set interval distance (every 50m, every 100m etc.)
Add from / to fields for chainages along the line.

I use ModelBuilder to automate this process as much as possible, but to generate the chainages, I need to use Excel because I have very limited experience with Python. So I generate the chainage values in Excel, then in the editor context in ArcGIS copy those values from Excel into the attribute table of my polyline features.

How can I transfer the logic I have in my Excel sheet over to Python to use in the Calculate Field tool within ModelBuilder?

The format of the from / to fields ends up being "K+MMM" where the K is the kilometer and the M's represent meters. So if I've got a starting point of 0km and 0m and an interval of 100m, the first chainage would read "0+000", followed by "0+100", "0+200".

Here are some screenshots of the logic in the Excel file: https://imgur.com/a/LHYyi0u


texts:
We have some polyline and point data in SHP format and would like to serve that data to an OpenLayers client.

First we thought about converting the data to GeoJSON and send it to the client as a whole package, but then we saw that the data is quite large (>500MB) and downloading such an amount as a whole is not a good idea.

So we came up with the idea to dynamically serve the data in a tiled format, using WMS.

We're not really experienced when it comes to do stuff like this, but at the same time, serving some SHPs as WMS shouldn't be too difficult.

So what is the easiest and smoothest way to do that? We've heard of GeoServer and MapServer, which one is easier to deploy (We're on CentOS 6). And is the caching of tiles included with serving WMS or does it have to be configured separately?


texts: How to make an editable City Map from OSM data?

texts:
I want to make a map of Cardiff, Wales, showing features useful to tourists, e.g. roads with road names, parks, museums, other POIs, etc. I want the final map to be clear to read - simple roads and names, basic outlines/boundaries - so editing (moving/manipulating) the data from its true georeferenced position would be ideal.

Using data from OpenStreetMap, I can't work out how to make a simple map in QGIS, and using Maperitive I can't find a way to edit the data from its original, simple view, despite being vector data.

A lot of the information online seems to be quite  old, for example the OSM plugin is no longer available in the latest version of QGIS.

Is anyone able to offer advise on making a simple, editable map of a city?


texts:
I am using OpenLayers 2.

I have a map with points from GeoJSON file.
I created a default style and a select style for these points.
I created a popup with attributes from GeoJSON when clicking on a point.
All works fine, except that I can't close the popup: the point is unselected, but the popup still remains even if I use the red cross.
I read this post that seems to deal with the same issue, but I can't use it in my application.

code for popup :

selection = new OpenLayers.Control.SelectFeature(operationCentroide);

    operationCentroide.events.on({
        "featureselected": monpopup,
        "featureunselected": onFeatureUnselect
    });

    map.addControl(selection);
    selection.activate();

    function onPopupClose(evt) {
        selection.unselectAll();
    }

    function monpopup(event) {
        var feature = event.feature;
        var contenu = feature.attributes.numope;
        popup = new OpenLayers.Popup.FramedCloud("chicken",
            feature.geometry.getBounds().getCenterLonLat(),
            new OpenLayers.Size(100,100),
            contenu,
            null, true, onPopupClose);
            map.addPopup(popup);
    }

    function onFeatureUnselect(event) {
        var feature = event.feature;
        if(feature.popup) {
            map.removePopup(feature.popup);
            feature.popup.destroy();
            delete feature.popup;
        }
    }


texts: How to get Legend Information from a WMS server?

texts:
I need to get the Legend Information from Geoserver in JSON or XML format, so that I can display it in my App in a custom way.

I'm basically looking for something equivalent to https://resources.arcgis.com/en/help/rest/apiref/index.html?mslegend.html in the OGC world.

The Image provided by GetLegendGraphics does not meet my requirements.

Is it possible to get Legend Information information in JSON/XML format?


texts: Adding Bing or Google imagery to the standard web application code derived from Geoserver

texts:
I am new to Geoserver. I know how to add layers to geoserver and view them using openlayers. Then i do a ctrl+U and copy the code into a html file which i place into my www folder in Geoserver. When i open this html file on the browser it works fine.
My problem is….
I want to add bing imagery or google imagery to this code, i am unable to make it work.
I have managed to use Bing imagery on a separate openlayers code, but, I am unable to integrate it with the code i copied from Geoserver. Do you have any tips on this?
I am not pasting my Geoserver code in this message because it is very long, and it is the standard code that one gets on pressing ctrl+U on the keyboard, when doing a layer preview in Geoserver.


texts: How to index and load a time-stacked NetCDF file with Open Data Cube?

texts:
Our use case is to make a time-stacked NetCDF file accessible via Open Data Cube.
Each NetCDF file (CMEMS data) has a time coordinate variable, with 8 time points evenly distributed over a range of one day.
When trying to datacube.Datacube().load() such a NetCDF file, the following error arises:
Error opening source dataset: NETCDF:"C:\CMEMS\mfwamglocep_2021013000_R20210122_mini.nc":VHM0
Traceback (most recent call last):
...
DeprecationWarning: Stacked netcdf without explicit time index is not supported anymore

How to achieve such an "explicit time index"?
Regarding indexing, it seems that such a time-stacked NetCDF file has to be covered by 8 time-specific dataset definition files (according to an issue discussion), but how to provide the specific time index?

Used resources:
• Open Data Cube version:
(cubeenv) C:\> datacube --version
Open Data Cube core, version 1.8.3

• Product definition:
name: mfwamglocep_mini
description: GLOBAL_ANALYSIS_FORECAST_WAV_001_027
metadata_type: eo

metadata:
  format:
    name: NETCDF
  instrument:
    name: NA
  platform:
    code: NA
  product_type: mfwamglocep_mini

storage:
  crs: EPSG:4326
  resolution:
          longitude: 0.08332825
          latitude: -0.08333588

measurements:
- name: VHM0
  dtype: float32
  units: m
  nodata: -32767
  add_offset: 0.0
  scale_factor: 0.01
- name: VMDR
  dtype: float32
  units: degree
  nodata: -32767
  add_offset: 180.0
  scale_factor: 0.01

• Dataset definition:
creation_dt: '2021-01-22T08:19:00Z'
extent:
  center_dt: '2021-01-30T13:30:00Z'
  coord:
    ll: {lon: -180, lat: -80}
    lr: {lon: 179.9167, lat: -80}
    ul: {lon: -180, lat: 90}
    ur: {lon: 179.9167, lat: 90}
  from_dt: '2021-01-30T03:00:00Z'
  to_dt: '2021-01-31T00:00:00Z'
format: {name: NETCDF}
grid_spatial:
  projection:
    geo_ref_points:
      ll: {x: -180, y: -80}
      lr: {x: 179.9167, y: -80}
      ul: {x: -180, y: 90}
      ur: {x: 179.9167, y: 90}
    spatial_reference: EPSG:4326
id: 272302c9-1449-4a33-8166-4b6083a8a711
image:
  bands:
    VHM0: {path: mfwamglocep_2021013000_R20210122_mini.nc, layer: VHM0}
    VMDR: {path: mfwamglocep_2021013000_R20210122_mini.nc, layer: VMDR}
instrument: {name: NA}
lineage:
  source_datasets: {}
product_type: mfwamglocep_mini
platform: {code: NA}

• Original NetCDF file, reduced to a publicly downloadable small subset (86 MB):
import xarray
ds = xarray.open_dataset("mfwamglocep_2021013000_R20210122.nc")
to_drop = []
for name in ds.data_vars:
  if name!="VHM0" and name!="VMDR": to_drop.append(name)
mini = ds.drop_vars(to_drop)
mini.to_netcdf("mfwamglocep_2021013000_R20210122_mini.nc")


texts:
Can there be some instances where it is necessary to interpolate DEM when re-projecting (e.g. SRTM from WGS84 to another projection)?

My current understanding of interpolation tells me that one should interpolate a dataset only once (i.e. as SRTM has already been interpolated there is no need to do that again).

I have come across some comments where apparently it has been interpolated (again) when re-projecting. This is just a general question, just want to know if my understanding of interpolation is right or wrong.


texts: Solving an Undefined Projection for US Counties (1990)

texts:
After inserting the cs99_d90.shp file obtained from the 1990 US Census, I received an error message claiming:

  The following data sources you added are missing spatial reference
  information. This data can be be drawn in ArcMap but cannot be
  projected.

After reading about several tools, such as the Define Projection tool which requires me to choose a projection from a list of projections that I cannot comprehend, I am hoping there is a much less convoluted solution to this problem. I have the 1990 County shapefile, and I want to layer it over a base map.

How can I solve this problem?


texts: Error in using/installing v.digit: GRASS GIS

texts:
When I wanted to edit a shapefile in GRASS I got this message:

  Unable to start wxGUI vector
  digitizer. Do you want to start TCL/TK
  digitizer (v.digit) instead?

  Details: No module named
  grass6_wxvdigit

Therefore, as suggested in this link:
http://osgeo-org.1803224.n2.nabble.com/wxpython-gui-and-v-digit-td3019439.html and
https://svn.osgeo.org/grass/grass/trunk/gui/wxpython/README - 6. VECTOR
DIGITIZER I have been trying to compile/install vdigit component. I could
install GRASS with C++ (--with-cxx), Python (--with-python) and wxWidgets
(--with-wxwidgets) support. But in the next step when I change to the
directory "cd gui/wxpython/vdigit" and run make I get the following error:

..............
cats.cpp:456: error: ‘Vect_rewrite_line’ was not declared in this scope

cats.cpp:461: error: ‘G_debug’ was not declared in this scope
cats.cpp:467: error: ‘Vect_destroy_line_struct’ was not declared in this scope
cats.cpp:468: error: ‘Vect_destroy_cats_struct’ was not declared in this scope

error: command 'gcc' failed with exit status 1
make: *** [OBJ.x86_64-unknown-linux-gnu/_grass6_wxvdigit.so] Error 1

I use Ubuntu 10.04 LTS with GRASS 6.4.1svn. I tried installing python-dev,
python2.6-dev packages as suggested in many forums but the error
continued to appear.


texts:
I'm struggling to find the math to do what I want.  I want to plot a top down view of a drillhole in ArcMAP 10.0.  I have the beginning UTM coordinate and then for every place in the hole that was measured, a dip (in negative degrees, horizontal being 0), a direction (in degrees, magnetic north being 0, although I can make this true north), and a length or depth.

I would just like to figure out the utm coordinates of each measured station so I can make a polyline.  This isn't math I know and although I'm finding lots of information about 3d vector math, I haven't been able to figure out how to translate it to my puproses.  So even if someone can tell me how to look for the formulas (what would it be called?) it would be much appreciated.  Once I know the formulas, I can use it in Python.  Thanks!


texts:
I have been given a DXF which is in mm instead of m. It also hasn't been spatially referenced.
I'm not too concerned with it being referenced in the right spot as I can easily translate it in another software package but is there anyway to convert it to m in ArcMap? I do have an advanced license.


texts: Editing X,Y of point in shapefile using ArcGIS for Desktop?

texts:
Using ArcMap 10.1, I edited the coordinates within a shapefile.

They are long/lat but the point will not move to that new coordinate on the map.

I've tried clicking on the Longitude column in the attribute table and 'calculating geometry' but still nothing occurs.


texts: Clipping raster with shapefile and exporting each raster subset using ArcPy?

texts:
I want to clip a raster image using a shapefile with differents polygons (each polygon have its "FID_malla" value [FID_malla is large data type]) and save the subset rasters individually with "FID_malla" name.

So if I have 50 polygons the result of the clip will be 50 raster images. To do that I use for iterator.

I run this code, but I obtain a error in SelectLayerByAttribute: Parameters are not valid  .
Also, I don't know if after SelectLayerByAttribute operation the code runs well.

def getParameterInfo(self):
    """Define parameter definitions"""

# First parameter
    param0 = arcpy.Parameter(
        displayName="Shape recorte",
        name="in_features",
        datatype="GPTableView",
        parameterType="Required",
        direction="Input")

# Second parameter
    param1 = arcpy.Parameter(
        displayName="Imagen satelite",
        name="raster",
        datatype="DERasterDataset",
        parameterType="Required",
        direction="Input")

    parameters = [param0, param1]
    return parameters

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True
def updateParameters(self, parameters):return
def updateMessages(self, parameters):return

def execute(self, parameters, messages):
    """The source code of the tool."""

a= parameters[0].valueAsText
b= parameters[1].valueAsText
#-----------------------------

with arcpy.da.SearchCursor(a, ["FID_malla"]) as cursor:
    for row in cursor:
        arcpy.SelectLayerByAttribute_management (a,  "NEW_SELECTION",  "FID_malla={}".format(row[0]))

        #Execute ExtractByMask
        out = ExtractByMask(b, a)

        out.save("C:\Users\mvinuales\Desktop\Image_Chips\imagen" + "\_" + str(row[0]))


texts: Comparison with different resolution rasters

texts:
The purpose is to compare the fractional snow cover of the coarser pixel to the binary snow cover (0 for no-snow and 1 for snow)  of lower resolution pixels. Both the datasets are of the same area and share the same WGS84 coordinates.
Below are two of the methods that I am thinking of:
As I understand random point sampling is independent of resolution in ArcGIS. So, I have two rasters one with coarser resolution say in kilometers and the other one is of finer resolution in meters. I create random points for each of the rasters. Then I run the "Extract multi values to point" tool to extract the pixel values of the attributes to these points. Please note that the attributes has to do with fractional area coverage in percentage (0-100%) in a pixel for one of the rasters, whereas the other raster which is of 500m resolution has binary/discrete values i.e. 0 /1 . The goal is to do a comparison between the point values based on fractional area coverage/percentage of each cell in the coarser raster to the binary data finer resolution raster.
So for example on a given day, the fractional covered area is 96% for a particular coarser pixel. As I know the coarser pixel has approximately 165 pixels, and so for that given day 60 of the finer resolution pixels have value 1. Hence, I add them and divide it by 165, basically ((60/165)*100)= 36.36%. I can then compare the two fractional area values i.e. 96% with 36.36%
So far I have found that there are 165 points representing the smaller pixels that lie in a single coarser cell. My question is somewhat similar to this question with the only difference that temporal resolution is same.
Will this be a valid comparison between the two sets of random point sampling, even though the resolution of the rasters is quite different?
Second Approach
I find the number of pixels in the coarser pixel, which is approximately 165, and then for those 165 pixels I add the binary one values of these pixels and divide them by total number of pixels which is 165, this might give me fractional area coverage of snow. I can then compare the fractional area coverage of the coarser pixel with the fractional area coverage of the 165 pixels. But even in this case, I am noticing that some of the pixels lie right on the boundary line of the coarser pixel.

The whole purpose is to comparison in such a way that one resampling of resolution can be avoided.


texts: Error when get the Z value to point from the attribute in QGIS

texts:
I have imported point from a csv file with the xyz coordinates in the attribute like this in QGIS.

And now I want to add Z value to the point. I try to use v.to.3d to add Z value but I get the error. How to fix this or Is there any tool to do like this?


texts: Leaflet WMS layers - Legend control not working

texts:
I'm having difficulties making my legend control work. I've loaded 3 WMS layers from Geoserver and declared their legends. I grouped the layers as Overlays and I have 2 basemaps. I found a code and adapted it as below to make a switch between the legends.What this does now is only displays Legend2.
In all the examples I've seen online, the layers (be it declared as baselayers or overlays) do not overlay in the app, but are removed upon another layer selection.

Does this function work only for layer removal or should it work in my case too, and I'm missing something?

var overlayMaps = {
  "DropOut_Rates": DropOUT,
  "UAT3":TAU_L3,
  "Schools": Schools
};

legend2.addTo(map);
currentLegend = legend2;

map.on('overlayadd', function (eventLayer) {

  if (eventLayer.name === 'DropOut_Rates') {
    map.removeControl(currentLegend);
    currentLegend = legend2;
    legend2.addTo(map);
  }
  else if (eventLayer.name === 'UAT3') {
    map.removeControl(currentLegend);
    currentLegend = legend3;
    legend3.addTo(map);
  }
  else if (eventLayer.name === 'Schools') {
    map.removeControl(currentLegend);
    currentLegend = legend1;
    legend1.addTo(map);
  }
});


texts:
I have raster layer imported in QGIS. How can I read value of the variable represented by the raster at mouse position? I just have a new raster layer and I want to read values at random places to have a basic idea about the raster.

I prefer not to create point layer for this purpose, but if it is needed, please tell me how to do it that way.


texts: Calculating Kernel Density within defined region?

texts:
I am attempting to calculate a Kernel density raster dataset from a point dataset within a defined regional boundary. However the output dataset generated is in a large square that cuts off one of the corners of this geographic boundary.

Is there a method to calculate the Kernel density within this polygon boundary?


texts: Spatial Join has Merge Rule Delimiter greyed out

texts:
I'm trying to use Spatial Join between a line layer and polygon layer. Where the polygon layer contains Area Names I would like to add to the lines. Sometimes a line will cross more polygons, so I'm thinking of using a Join_One_To_One, with the Merge Rule as Join for names. This works fine however, I would like a delimiter in the Merge Rule, such as a Comma.

However the Delimiter part is greyed out, any suggestions why and how to get around it?


texts:
I am writing below SQL query in mapbasic to insert columns from two tables t1 and t2 into table Sites but error is occuring "Variable or field t2.CELLID not defined"

Table Sites (Source integer, N integer, LAT float, LONG float)
Select t1.CELLID, t2.CELLID, t2.LATITUDE, t2.LONGITUDE from t1 left join t2 on t1.CELLID=t2.CELLID where t1.obj within zone into Sites

Can anyone tell me how to rectify it?


texts: Converting Microsoft Visio 2007 .vsd file to xml workspace document (without the macro)

texts:
I was just editing my geodatabase schema in Microsoft Visio and noticed that it is the new 2007 version, rather than 2003 which had a handy export tool, therefore I downloaded the ESRIExportToXML macro, however when I use it, it comes up with the following error:

Run-time error'-2032465751 (86db08a9)':

The line of code which it breaks on is the following:

Set pAO = Application.Addons.Item("UML Background Add-on")

I have tried to fix the code or find other tools to export it to xml from visio but have virtually given up.  Therefore I was wondering if I could do something along the lines of export it to pdf then, using acrobat pro save as an xml (I have tried this but it created an invalid xml file), or convert it in excel?  Therefore is there a way round or alternative route for converting an .vsd file to xml file?


texts: While clipping a Raster generated from SCP plugin (classification) all the colors are vanished in the end output and the raster generated is in B/W

texts: I have done all the LULC classification using SCP plugin. after completing my classification, I wanted to clip the raster using a Shape File. But when I do it, the colors are vanished in the output clipped raster. The raster generated is in Black and White. Help me to fix this.


texts: Creating polyline templates (polyline vs polyline ZM) for future network dataset?

texts:
I am about to request large amounts of polylines from a variety of sources (municipalities). The idea is to collect bike path data and build a regional bike network (I work for the region administrative board).

To help municipalities (and myself once I need to assemble all the data), I plan to create template files. I now realize I do not know enough about network datasets and 2D to 3D conversions. I basically don't know what I don't know, hence this very open question.

I have:

Arcgis 9.3 (hopefully 10.2.1 in the fall but not for certain).
High resolution DEM of the region.
limited bike path shapefiles (segmented polyline ZM, not a network
dataset)

I need:

advice in creating template shapefiles that are easy and flexible enough for the municipalities to use and still be fairly straight forward for me to put together for a network.
The templates will contain 3 attributes (short integers) that denote road type, traffic situation and recreational value.

specifically I would like to know more about polyline vs polyline ZM? - bear in mind that the municipalities are using a variety of software (CArcGIS, MapInfo, CAD(?) etc. etc.).


texts:
Is there any information out there on how to start creating scripts using Python and QGIS 2.10? I only have a little coding experience and an very unsure of how to start.

My first script would probably be something like a process for creating a buffer from a point, setting the properties of the buffer too if possible (e.g. orange, 0.8 line). I do simple tasks in QGIS often and I think scripts would help me a lot.


texts: How do I iterate a clipping process through a list of geodataframes?

texts:
I have a folder full of shapefiles, I would like to read them into geopandas and then clip them all to one boundary iteratively. I have loaded all of the shapefiles into a list with this correctly thanks to @Kadir Şahbaz:
import glob
import geopandas as gpd
shapefiles = glob.iglob('E:/folder/shapefiles/*.shp')

gdfs = [gpd.read_file(file) for file in shapefiles]

However now I am stuck on iterating through each shapefile (geodataframe in the list gdfs). The code below only applies the clip to the first shapefile(item) in the list:
for shapefiles in gdfs:
    clipped = [gpd.clip(shapefiles, boundary)

I need to iterate the clip through each shapefile(geodataframe) in the list (and later save each as a separate shapefile). Any help would be greatly appreciated!
Thanks


texts: Formatting street network data for ArcGIS Network Analyst?

texts:
How do I create a good street network for network analysis?

How should I format the street data to import it into Network Analyst, including prefixes, directions, intersections, and geocoding addresses?


texts: What is ID in leaflet code for L.tileLayer

texts:
I'm working with this famous code:

L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}',
    maxZoom: 18,
    id: 'CAN'T FIND IT',
    accessToken: 'FOUND IT'
}).addTo(mymap);

I can't find the ID of my mapbox but when i change it to mapbox.streets, the map shows up.
Mapbox gives me a style URL, what should i do with it?


texts:
I'm looking for some application to view some large TIFF images and large mosaics consist of tiff images.

At the moment I have VLIV for viewing large tiff files.
And tools based on HD view: HD View Utilities and Image Composite Editor.

Any suggestions from GIS world?


texts: Rule-based symbols and Expression String Builder

texts:
I am working on a data set which is looking at total number of fish landings in accordance with vessel length within the North Sea. After plotting the data I want to be able to categorize total number of fish landings in accordance to the length of the vessel.

I have attempted to use the Expression String Builder within the Rule-based symbol in properties and symbology, in which there are four categories available <=10,>10-<=12,>12-<15, and >=15. I want to be able to plot the length LENGTH_GP of the vessel to total number of fish landings, with each category displaying a different colour.

Does anyone have any suggestions to this issue?


texts: How to export a .sid image as high DPI .tiff or .bmp from QGIS?

texts:
I have 2 50 cm resolution areal photographs in .sid format.
How can I export these layers as a regular image format (like tiff or bmp) which I can then open in any simple picture viewer program?

I've tried saving a map, but that reduces the image resolution. I need to resulting image to retain it's original resolution.


texts: Data Source: Address ranges for every block in NYC

texts:
I'm looking for data that defines the address ranges used on NYC blocks. On this screenshoot you can see google has noted the first addresses on each block of West End Avenue and 92nd/93rd/94th st.  (658, 680 and 700 respectively).  Based on that info I can extrapolate assuming that there's an even/odd side of the street and there's no unrepresented addresses to get something like this (inefficient representation used only as an example):

Street | btwn1 | btwn2 | first | last | side
-------------------------------------------------
 W End | w 92  | w 93  |  658  | 678  | west/even
 W End | w 92  | w 93  |  659  | 679  | east/odd
 W End | w 93  | w 94  |  680  | 698  | west/even
 W End | w 93  | w 94  |  681  | 699  | east/odd

Can anyone suggest a source of this type of data? Basically it's a poor mans geocode so I can convert an address (Jeryl Metz Books - 697 West End Ave) into a suggestion of how to get there (West End btwn 93/94) even without an internet connection.  I'm sure this data is contained within OpenStreetMap, but I don't have a clue as to how one might extract/generate what I need.

.


texts: Batch saving MXDs to group layer files using ArcPy?

texts:
Does anyone have / know of a script that converts all the MXD`s in a folder to individual group layer files ?

It needs to use the MXD name as the group layer file name.


texts:
I did a module in grass wich I called r.out.haralickFeatures.py. I did a .qgm. I can't get the parameters rigth.

#!/usr/bin/env python

#%module
#% description: Compute asm, contrast, var, idm and entr for each band of a TIFF file
#% keyword: raster
#% keyword: haralick
#% keyword: features
#%end

#%option G_OPT_F_BIN_INPUT
#% key: input
#% type: string
#% required : yes
#%end

#%option
#% key: size
#% type: integer
#% required: yes
#% multiple: no
#% key_desc: value
#% description: The size of moving window (odd and >= 3)
#% answer: 3
#%end

#%option
#% key: distance
#% type: integer
#% required: yes
#% multiple: no
#% key_desc: value
#% label: The distance between two samples (>= 1)
#% description: The distance must be smaller than the size of the moving window
#% answer: 1
#%end

#%option
#% key: categories
#% type: integer
#% description: Number of categories to rescale/recode the image
#% answer: 255
#% required : yes
#%end

#%flag
#% key: r
#% description: Rescales instead of recoding the image
#%end

#%option G_OPT_F_OUTPUT
#% key: output
#% type: string
#% required : no
#%end

import sys
import grass.script as gscript
from grass.pygrass.modules import Module, ParallelModuleQueue
from copy import deepcopy
from osgeo import gdal
import os
import time
import atexit

def main():
    options,flags = gscript.parser()
    inputFile = options['input']
    sizeMovingWindow = options['size']
    distance = options['distance']
    numCategories = options['categories']
    output = options['output']
    recode = flags['r']

    print options

    ...

if __name__ == '__main__':
    atexit.register(cleanup)
    main()

and the .qgm:

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE qgisgrassmodule SYSTEM "http://mrcc.com/qgisgrassmodule.dtd">
<qgisgrassmodule label="Compute asm, contrast, var, idm and entr for each band of a TIFF file" module="r.out.haralickFeatures.py">
    <file key="input" />
    <option key="size" answer="3" type="integer" />
    <option key="distance" answer="1" type="integer" />
    <option key="categories" answer="255" type="integer" />
    <flag key="r" answer="off" />
    <file key="output" type="new" filters="Geo TIFF (*.tif *.tiff)"/>
</qgisgrassmodule>

The output i get after running it in qgis is:

C:/OSGEO4~1/apps/qgis-ltr/./grass/scripts/r.out.haralickFeatures.py
input=C:/Users/Carolina/Documents/PPS/Informe/Landsat8-OLI.tif size=3 distance=1 categories=255 output=
{'input': '@ARGS_PARSED@', 'output': '', 'distance': '1', 'categories': '255', 'size': '3'}

ERROR: Unable to open file:  @ARGS_PARSED@

Successfully finished

Anybody knows why i get @ARGS_PARSED@ and not the parameters i'm supposed to get?


texts:
My version : postgres 8.4.1/postgis 1.4.0

I want to perform a simple spatial join with the following query :

SELECT code_insee, nom, id_canton, lib_canton
FROM localisants.alca_toponyme_commune
LEFT JOIN localisants.alca_pseudo_cantons
ON ST_Intersects(alca_pseudo_cantons.shape,alca_toponyme_commune.shape);

The type of geometry is POINT for alca_toponyme_commune and
MULTIPOLYGON for alca_pseudo_cantons.
Both tables are planned in EPSG:2154.

In my result, the join fields (id_canton, lib_canton) are empty.

Is there a problem in my query?


texts: QGIS: Cannot add negative numbers in a field

texts:
I am using QGIS 3.4 and when I finish drawing a feature (line), I try to add a negative value (-1) in the attribute dialog and it is automatically changed to NULL. The field is integer.

Even if I leave it NULL and then open the attribute table and try to edit it there, it still does not let me use a negative number - automatically changes to NULL. The only way is to change it using Field Calculator and update the existing field to -1.

I tried the same in version 2.18.17 and it works fine.

I haven't used QGIS for years and I am not entirely familiar with the changes in the latest versions. Is there a new setting/option that I don't know of?

EDIT :

I have reported the bug and it has now been fixed.

PostGIS/geojson: Cannot input negative numbers in int fields

The new QGIS versions won't have this issue


texts: Removing contour artifacts from low slope regions

texts:
I am generating contour lines from a DEM using QGIS and the python console. The process generates many small circles in regions of low slope (little black circles) that detract from the long connected contour representing the "true" contour (indicated by hand in blue) in the picture.

As suggested "https://community.esri.com/thread/206552-isyhypses-filtering-out-very-small-areas-of-elevation" I could filter the contour lines for length but that may also remove many "true" contours that can be seen in my attached image such as the few small hills here and there (indicated by hand for example in red).

Google searches suggest to smooth the underlying source DEM as in smoothing DEM but that would reduce the quality of the generated contour lines.

The best idea I had, would be to somehow boolean filter the contour lines against length as described above AND/OR local slope % as described calculating-slope-in-percentage-with-qgis, but I haven't yet determined how or if this is possible in pyQGIS - possibly by calculating a slope % tif and then reading it at every vertex of each contour then averaging for each contour feature to create a second attribute to boolean filter along with length. Before I embark on that (probably time consuming)  solution, is there an existing method for removing the low-slope artifacts generated in the DEM to contour process available in pyQGIS?

Edit: image added to indicate results of resampling (1s blue, 18s green, 27s red, underlying DEM 1 arc second)


texts: Extracting position date of max and min value in NDVI time series in Google Earth Engine?

texts:
I need to extract the DATE for max and min NDVI values in a MODIS NDVI time series.
i think I am able to extract the date for max but I have problems with the minimum.

// turn image collection into an array
var array = NDVI.toArray()

// sort array by the first band, keeping other bands
var axes = { image:0, band:1 }
var sort = array.arraySlice(axes.band, 0, 1);  // select bands from index 0 to 1 (NDVI)
var sorted = array.arraySort(sort);

// take the first image only (MAX NDVI)
var length = sorted.arrayLength(axes.image)
var values = sorted.arraySlice(axes.image, length.subtract(1), length);

// convert back to an image
var max = values.arrayProject([axes.band]).arrayFlatten([['ndvi', 'time']])

// get the max band
var ndviMax = max.select(0) // ndvi max
var time = max.select(1) // get day

This is the link to the code
https://code.earthengine.google.com/29618e626ddc93d7d6a6d00103518c6e

Any help/suggestion?


texts: I'm wondering if there's some way to find the orientation of a set of points along a line as a number (degrees from north) in ArcGIS?


texts: Location Data in TAIP Format to Geoevent Processor

texts:
We are testing whether we can get the surveyed cordinates in a Trimble GPS transferred online via TAIP format ( Trimble exchange protocol) to ArcGIS Geoevent processor .

For this , one main concern is how the device to configure the Device to receive the data in TAIP format and secondly how to push the TAIP format to the Server.

The page below provides some insight to the work flow.

http://www.arcgis.com/home/item.html?id=6b458b1f4de5436486ffbd7801c6b438

The question is , how can I override the programs in the GPS to configure the TAIP format . I guess this qn should be put to the GPS hardware dealer, but if anyone can offer a suggestion, we make some headway..


texts: Fetching the leaflet map with plain JSON data

texts:
countnuing the problem from this query:

Import pure (regular) JSON file to Leaflet

I have trouble with input the plain .json data to the leaflet map. I tried 2 ways to do this:

I attached the .json code as a JavaScript variation:

  var bv = {
 "Sheet1": {
"-430524279882": {
 "MDU_ID": -430524279882,
  "Status": "Invoiced",
 "THP": 12,
 "Latitude": 52.41635,
 "Longitude": -1.55227
  },
   "-430532279893": {
  "MDU_ID": -430532279893,
  "Status": "Invoiced",
  "THP": 12,
  "Latitude": 52.41635,
  "Longitude": -1.55227
   },
   "-430622279798": {
   "MDU_ID": -430622279798,
   "Status": "Invoiced",
   "THP": 12,
   "Latitude": 52.41555,
  "Longitude": -1.55155
   },
    "-430519279894": {
   "MDU_ID": -430519279894,
    "Status": "Invoiced",
   "THP": 12,
   "Latitude": 52.41635,
   "Longitude": -1.55227
    },
    "-430558279922": {
   "MDU_ID": -430558279922,
   "Status": "DO NOT DESIGN - SDN",
   "THP": 5,
   "Latitude": 52.41635,
   "Longitude": -1.55227
   },
   "100071574860": {
   "MDU_ID": 100071574860,
   "Status": "DO NOT DESIGN - SDN",
   "THP": 5,
   "Latitude": 52.41555,
   "Longitude": -1.55155
    },
    "-430564279884": {
    "MDU_ID": -430564279884,
    "Status": "DO NOT DESIGN - SDN",
    "THP": 5,
    "Latitude": 52.41635,
    "Longitude": -1.55227
   },
    "-435086278038": {
    "MDU_ID": -435086278038,
   "Status": "Invoiced",
   "THP": 10,
   "Latitude": 52.39908,
    "Longitude": -1.4861
    },
    "-476625255228": {
    "MDU_ID": -476625255228,
    "Status": "Invoiced",
    "THP": 6,
    "Latitude": 54.18998,
    "Longitude": -3.88052
   }
 }
};

and next:

   var myPointsLayer = L.featureGroup(bv);
   for(var key in data2.Sheet1){
   var latLng = L.latLng([data2.Sheet1[key].Latitude,
   data2.Sheet1[key].Longitude]);
   L.marker(latLng).bindPopup(data2.Sheet1[key].Status).addTo(myPointsLayer);
  }

My 2nd way was:

var url = 'bv.json';

var laya = $.getJSON(url);
var myPointsLayer = L.featureGroup(laya);
for(var key in data2.Sheet1){
  var latLng = L.latLng([data2.Sheet1[key].Latitude, data2.Sheet1[key].Longitude]);
  L.marker(latLng).bindPopup(data2.Sheet1[key].Status).addTo(myPointsLayer);
}

Both not working:

Console says:

> Uncaught ReferenceError: data2 is not defined
>     at (index):119

whereas data2 belongs to my .json file.

Is there any way to populate this stuff on the map?