texts: Remove internal boundaries in spatial object/polygon for plotting

texts:
I am trying to plot a spatial object where only the outer boundaries/borders are plotted. The only way I can think to do this is by manually removing all the lines/points within the spatial object and leaving those for the outer borders in it, but I have no idea how I know which are which and wondered if there was a simpler way to do this.

The file I am using can be found here

The code I use to read this in and plot is:

m1<-geojson_read("topo_E06000023.json", what="sp")
plot(m1)

The object plotted is shown below and currently this shows all the lines/borders.


texts: SpatiaLite : update table with a spatial relationships

texts:
In a SpatiaLite DB, I am trying to update the value of a set of entities of a layer (table1) selected by a spatial relationship (equals) with entities from another layer (table2)

I can perform a simple selection that take about 1 minute with this request :

SELECT table1.* FROM table1, table2
WHERE EQUALS(table1.geom, table2.geom);

I can perform the update but it takes about 1 hour with this request :

UPDATE table1
SET field1 = 1
WHERE table.id IN (
SELECT table.id FROM table1, table2
WHERE EQUALS(table1.geom, table2.geom));

When I try this request which I hope could take the same time as the first request :

UPDATE table1
SET field1 = 1
WHERE EQUALS(table1.geom, table2.geom);

I get this error :

SQL error: no such column: table2.geom

I wonder if I can do the update faster than the second request. If I do it with QGIS, first the selection, then the update it very fast but I need to do it with SpatiaLite.


texts: Location of UK Town to Postal Sector database or .csv

texts:
I've been searching all over, but cannot find a database/CSV for matching a UK Town to postal sectors.

Example:

Hampton Hill” is postcode sector” TW12 1


texts: Importing information on bridges into postgis-table using OSM2PGROUTING

texts:
In OpenStreetMap (OSM) road data there is information about whether or not a road feature is a bridge or not (boolean). However, when importing into a postgreSQL data base for routing (pgrouting), using osm2pgrouting, I cannot find information whether a feature is a bridge or not. Is it possible to tell osm2pgrouting, or osm2po, to include this information somehow? Maybe in the mapconfig.xml?

I understand that OSM2PGROUTING uses the bridge-information for excluding nodes at a bridge-road intersection, but I need that information in later analysis.

Here is the column-heads for OSM road data:

Versions:

Windows 10
Postgis 2.2.1
postgreSQL 9.5
pgrouting: 2.1.0
OSM2PGROUTING 2.1.0


texts:
I’m fairly new to QGIS.
I have a polygon that I need to split up into circles. I have the following requirements:

I must be able to specify the radius for these circles, for example 2 km, so that they’re all the same size.

The circles may overlap, but should do so as little as possible.

Together the circles must cover the entire polygon, and it’s preferred that they cover as little of the area outside the polygon as possible.

Is there a way of doing this?


texts:
I'm looking for best practices (plugins, procedures...) in QGIS to generate an HTML document (which later on can be finalized in MS Word) in which the result of spatial queries can be documented.

An example: I have a layer with parcels of a municipality. I would like to generate a document, according to a pre-defined HTML layout, where the results of overlays of this parcel with other geographic layers (e.g. legal boundaries) are included. E.g. in my document, there is a question "is the parcel located within an agricultural area?" and behind the question there are two checkboxes (yes and no) which should be checked depending on the result of a spatial overlay between the selected parcel and the layer with agricultural areas.

In addition to the above-mentioned spatial query, it would be also useful to include in the HTML-document the information from some attribute fields from the layer with agricultural areas (e.g. name of the agricultural area) in case the parcel overlaps with an agricultural area?

Are there any easy-to-use plugins or tools for this in QGIS?


texts: Find parcels dupliactes in different street names

texts:
I want to find duplicates in parcels, deep in streets basis.

My schema goes like this; Neighborhood >> Block >> Parcel >> Street

Let's say I have Street1, Street2 and Street3, all of them might have duplicated parcels. I want to do it through a model, to search all the streets "individually" then give me results of duplicates if there any.


texts: Evenly spaced labels for polyline layers

texts:
I'm using QGIS 3.8 and am attempting to label roadway segments (OSM) of varying length. By checking the "Merge connected lines" box and hiding features below a certain size, I can generally limit the display of labels.

However, because of sharp turns in lines, or places where there are two parallel lines for the same segment, overlapping or doubling still occurs. Furthermore, in areas where all segments are under my specified size, there's no label at all.

I want to label all features, but no two labels of the same value (e.g. "35") should be within some distance of each other. Is that possible?


texts: Storing information for each segment of spatial line

texts:
I have some spatial lines from a shapefile representing roads. I'd like to attribute to each segment an information, for example the state of the segment or more generally a scalar that gives an information on the segment. For example I'd like to map portion of the roads that need to be repaired or the portion of the roads that no longer exist.

As far I know, we can only attribute values to individual lines and this prevent to map anything that applies only to some portion of the lines.
My question is language and software agnostic. I'm actually able to deal with that programmatically using lines and points but them I don't know how to share that in a simple way that is compatible with GIS software.
Is there some file formats supported by spatial software to achieve this task?


texts: Split Polylines and Count Number of Segments Tool in ArcGIS Pro?

texts:
I recently stumbled upon a tutorial on Esri's website that would split up a polygon or series of polylines (can't remember specifically) but I lost the link!
Is anyone familiar with a tool like this?
There was a very specific diagram that I tried to recreate in photoshop:

The result of the operation would be something along the lines of 8 lines being counted from this one shape.


texts: Creating polyline-based "heatmap" from GPS tracks?

texts:
This winter I am planning to track my downhill skiing/snowboarding using a GPS.  Most of my riding will occur at the same resort.  I would like to be able to create a sort of "heatmap" that shows the amount of runs that I have made in a given area.  As I add more and more GPS traces to my database, my goal would be to see a sort of linear heatmap of the most traveled areas.  Given the nature of downhill skiing, you would expect the uphill chairlift lines will be the "hottest" because it will be the only places visited over-and-over again.

Given that 1) my track will not be the same every time and 2) the area covered by following one "run" may be a few hundred feet wide, what might be the best way to analyze this "linear" data to create a sort of heatmap?  My thoughts were to buffer the lines, then intersect the polys to get a sort of Venn Diagram thing going.  My preference is to use open-source technologies.  I've got QGIS and PostGIS loaded and available.

UPDATE: In regards to @blah238's response, I was thinking of something that might be able to "collect" the number of passes ("runs") through an area, and then symbolize by the count.  Conceptually, this would be similar to ArcGIS "Collect Events" (but for lines, not points) or Collapse Dual Lines To Centerline (but for multiple line in roughly the same area).

A more visual example of a similar concept might be a traffic-flow map, where highly-congested areas would equate to "highly-traveled" ski runs/areas:

I've read the following questions that may give some ideas, but they don't really address what I am trying to accomplish:

Clustering Trajectories (GPS data of (x,y) points) and Mining the data

Managing error with GPS routes (theoretical framework?)


texts: Is there any way to display a number as currency in the hover window? I'd like the  value to read “$5.6 million” rather than “5555555,” for example.


texts: Reading an image as array in python using GDAL?

texts:
img=gdal.Open("D:\data\sub_66")
inputArray=img.ReadAsArray()

gives Error:
ValueError: array is too big.
What is the problem with this?


texts: Creating list of feature layers, running select by location, and returning results in ArcPy?

texts:
I want to perform Select by Location with a number of point layers and a corresponding number of polygon layers. For example, if I have a list of point layers (pt1, pt2, pt3) and a list of polygon layers (pg1, pg2, pg3), I want to know how many points in pt1 intersect pg1, and the same for pt2 and pg2, and pt3 and pg3. I eventually want to be able to turn this into a tool that can be used on various datasets.

I have been working on a Python script to accomplish this, and the first half works great, but the second half, not so much.

#Import arcpy and set working environments
import arcpy
arcpy.env.workspace = r"Z:\Scratch\ChristinaScratch\pythonpractice.gdb"
arcpy.env.overwriteOutput = True

#Create list of only point feature classes that are present within the
workspace and convert to feature layers

pointFeatureClasses= arcpy.ListFeatureClasses("*","Point", " ")
for fc in pointFeatureClasses:
    arcpy.MakeFeatureLayer_management(fc, fc + "_pfl", "", "","")
    print fc
areaFeatureClasses= arcpy.ListFeatureClasses("*","Polygon", " ")
for fc in areaFeatureClasses:
    arcpy.MakeFeatureLayer_management(fc, fc + "_afl", "", "","")
    print fc

#for loop to iterate through each point and each area feature class in the
list, select the points that intersect with the area feature class,
and return the number of records selected

pointFeatureLayers = arcpy.ListFiles("*_pfl.lyr")
areaFeatureLayers = arcpy.ListFiles("*_afl.lyr")

for fl in pointFeatureLayers:
    print fl
for fl in areaFeatureLayers:
    print fl

for point in pointFeatureLayers:
    for area in areaFeatureLayers:
        arcpy.SelectLayerByLocation_management("point","INTERSECT","area")
        result = str(arcpy.GetCount_management(point))
        count = int(result)
        print count

I am able to create the list of feature classes and convert them to feature layers without problem, but the rest of the script doesn't seem to be doing anything. I'm not thrown any errors, I'm just not getting anything. I just discovered that while the created feature layers are added to the map and are sourced back to my workspace gdb, when I look into my gdb through ArcCatalog they aren't actually there.

What is going wrong?

I'm fairly new to Python.


texts: How to show paths in Google Earth on mouse-over?

texts:
I've seen that paths included in a kml file can be configured to show up only when you do a mouse-over:
http://www.barnabu.co.uk/flight-route-maps-for-australia-and-new-zealand/

How can that be done?


texts:
I am working with a DEM that has a bit more data than I actually need.  Specifically, even the lowest resolution version is 100ft and I'm trying to reduce it down to 500ft or even 1000ft intervals.  I should disclose here that I am pretty much a total noob to both ArcMap and geodata generally at this point.

After digging around on the internet, I think the steps to doing this are:

Double click on the 100 Foot layer in the table of contents window
Click on the definition query tab
Click the Query Builder button
Create a query like ELEVATION%500 = 0
Profit

Unfortunately that doesn't seem to work.  Taking a look at the data (right click on 100 Foot layer -> open attribute table) it appears that the elevation for each layer is not simply a number ("0").  Instead, it is a number and the unit of measure that every number in that column uses ("0 ft").  I assume  this is what is preventing me from coming up with any matches.

If I am, in fact, on the right path here, is there a way to construct a query that ignores the " ft" portion of the data (or just strips out " ft" from the field)?  Alternatively, if I am not, in fact, on the right path here, does anyone have any suggestions for a potentially more fruitful line of inquiry?


texts: Getting 001270 error : Consolidating the data failed while publishing as map service in ArcGIS for Server?

texts:
I am getting an error "001270  error" .

I am inserting table for sql-server database and joining few records to polygon layer

also I have lat,Lon data in few tables so I am displaying data using "Display x,y values" and publishing all data into arcgis server 10.1 (using ArcGIs Desktop 10.1)

While publishing I was getting above error..

also my database & folders are registered in ArcGIS Server (through ArcCatalog I did this)

I found following article which explain that its bug  Link : http://support.esri.com/en/bugs/nimbus/TklNMDgzMjc5 ..but I am not sure.


texts: QGIS: zoom further than threshold in WMTS

texts:
Is there a way to zoom further than the threshold allowed by a wmts feed?
Of course the idea is not to download any more images as that would be impossible. I just want to zoom on the already downloaded images.

I'm using IGN's BDORTHO with the "pratique" key and as soon as I zoom further than some threshold (which strangely varies depending on the computer I'm using it on), the map goes blank instead of scaling the image.


texts: How to remove nodata values from XYZ output with GDAL?

texts:
I have a buffered line route that is overlayed on a DEM:

I've clipped the DEM to the buffered line, and exported it as an ADF ESRI grid file. I then used gdal_translate to convert it from the Arc/Info ASCII Grid to the ASCII Gridded XYZ. This all works fine, but the output XYZ has tens of thousands of entries for elevation marked as -32768 which I assume is nodata.

How can I have the output XYZ file only contain values that aren't -32768? Is there a tool or parameter in GDAL I'm missing?

And why are there nodata values? Is the raster I export after clipping really a rectangle with nodata values where the buffered line isn't present?

The XYZ file is going to be used in PLS-CADD.


texts: Working around US Census site blocking of IP addresses?

texts:
I'm downloading TIGER data from the census site for the TIGER geocoder,  but the computers IP has been blocked.

I'm trying the following now:

Generate the loader script on the target machine.
Copy that file into a different machine that is not banned.
Keep only the wget statements in order to download the data.
copy the downloaded data to the target machine.
Edit the loader script commenting wget statements.

Is this a feasible approach?


texts: Adding incremental id value to features with certain attribute value

texts: I am using QGIS 3.12. I have a layer of trees and a virtual field 'tree number' that counts all trees that I add to the layer ($id). I want to set that if a height of a tree is more than 3 (in field 'height') it will not be counted in the 'tree number' field. Can it be done?


texts: Hole in tiles with GeoServer 2.18 (while it works in 2.13)

texts:
I'm upgrading my system from GeoServer 2.13.0 to GeoServer 2.18.2 on docker using kartoza image.
All works fine just I have this issue with a raster layer based on several pyramidal rasters: I get holes on some tiles, not on all of them, and sometimes only on one color band, like in this picture:


texts: How to develop responsive layer container in Leaflet?

texts:
I'm using leaflet.js to develop a drawing layer tool with the full-page map. Users will draw a group of layers to describe an industrial product's parts.

The leaflet map and its layers should be fit to screen. For example, if a user draws a circle to the left bottom corner with 1920x1080 resolution pc, I should see it in the left bottom corner with my 1366x768 resolution pc. The leaflet map is resizing with the page correctly but the circle is being out of the screen in my pc because of its latitude-longitude in this example.

I tried to develop this example with invalidateSize, maxBounds, fitBounds, etc. but haven't succeeded yet. How can I develop this scenario?

Map (React):

  const corner1 = L.latLng(-180, -180);
  const corner2 = L.latLng(180, 180);
  const bounds = L.latLngBounds(corner1, corner2);

        <Map
        ref={map => {
          this.tryToCreate(map);
        }}
        className="c-leaflet-map-map"
        center={[0, 0]}
        zoom={0}
        maxZoom={10}
        zoomControl={false}
        onzoomend={this.zoomLevelsChange}
        maxBounds={bounds}
        maxBoundsViscosity={1.0}
        zoomSnap={0}
        crs={L.CRS.Simple}
      >
        <ImageOverlay
          bounds={map.backgroundImageBounds || map.bounds}
          attribution="&amp;copy TofnaTech"
          url={
            map.backgroundImage ||
            map.backgroundImageURL ||
            TransparentImage
          }
          opacity={map.opacity}
          className={cn({ transparentImage: !map.backgroundImage })}
        />
        <FeatureGroup
          ref={controller => {
            this.controller = controller;
            this.tryToCreate();
          }}
        >
          <EditControl
            position="topleft"
            onCreated={this._onCreated}
            onEditStart={this.onEditStart}
            onEditStop={this.onEditStop}
            onEdited={this.onEdited}
            draw={{
              polyline: false,
              polygon: {
                shapeOptions: { ...defaultShapeOptions }
              },
              rectangle: {
                shapeOptions: { ...defaultShapeOptions }
              },
              circle: {
                shapeOptions: { ...defaultShapeOptions }
              },
              marker: false,
              circlemarker: false
            }}
            edit={{
              remove: false
            }}
          /></Map>


texts: Getting length data from GPX files within QGIS

texts:
I just got a Garmin etrex and promptly took it for a hike.  I can get the path and the trackpoints into QGIS with no issue.  However, I cannot figure out how to get the distance between the trackpoints.  I've played with the CRS and the file formats (converting to .shp for example).

When I go into the attribute table and attempt to calculate a new field using $length, all I get are NULL values.

Any ideas?


texts: I'm trying to parse up the addresses in a table using the above code. It runs but only produces NULL. I have run out of ideas


texts: I want to select counties along a river. The following picture is an output image by using the Ohio River. The river vector is a line or a buffered area. Now, I choose "Edit" and "Select Features" by hand. I would use "Select by Expression." But I don't know how to write "Expression" in two layers.


texts: Using shp2img and Remote Objects instead of WMS

texts:
I am a Flex developer that works with ESRI's REST Flex API most of the time, but have been exploring other options in the Flex world of web mapping.  I have been using WebOrb's implementation of Remote Objects for other, non-spatial development and have been totally blown away by the performance gains I have seen compared with traditional web services (SOAP and even HTTP).  And this has given birth to the following idea.

Instead of setting up a traditional WMS with MapServer, I am toying with the idea of firing the shp2img directly via remote objects, and sending that image back to the client to be drawn on top of the map. The only drawback to this approach that I see is that it is not very interoperable, i.e., no one else can consume this data but I am not too concerned with that right now.

I am comfortable with MapServer but far from a power user, so I am hoping to get some potentially valuable insights(i.e. warnings,ridicule...etc) from other MapServer experts.  This idea is based on the assumption that shp2img will be faster at generating a tile than its equivalent WMS GetMap request.

Any thoughts.


texts: Creating a heatmap from vector grid and point layer with QGIS

texts:
I would like to make a heatmap based on the number of hits within a certain area.

I begin by importing the data of my hits:
Layers > add layer > add delimited text layer = ok
I add google satellite layer = ok
I add a raster:
Vector > analysis tools > vector-grid > Select area from map > X=10, Y=10 > grid as polygon
Vector > analysis tools > count points in polygon
Open attribute tabel = I get my dimensions X min/max and Y min/max + number of points

Now I would like to make a heatmap where each square gets separate color. I think I need to use the data from the atttribute table but I don't know how.
I would like to obtain the result below

Example: each square with hits between 10-15 hits = red, 5-9 hits = orange, 1-4 hits = green.

Can this be done and if yes, how?


texts: Using GeoPandas intersections to fit a grid on a geometry

texts:
I'm trying to write a function that creates a grid over a vector shapefile using GeoPandas:
def grid(src_vector, out_vector, size):

    # read the vector data in mercator proj
    gdf = gpd.read_file(src_vector).to_crs('EPSG:3857')

    # extract bounds from gdf
    min_lon, min_lat, max_lon, max_lat = gdf.total_bounds

    # compute the longitudes and latitudes top left corner coordinates
    longitudes = np.arange(min_lon, max_lon, size)
    latitudes = np.arange(min_lat, max_lat, size)

    # create the grid centroids
    points = []
    for coords in product(longitudes, latitudes):
        points.append(Point(coords[0], coords[1]))

    # create a buffer grid in lat-long
    grid = gpd.GeoDataFrame({'geometry':points}, crs='EPSG:3857') \
        .buffer(size, cap_style=3) \
        .intersection(gdf) \
        .to_crs('EPSG:4326')

    # filter empty geometries
    grid = grid[np.invert(grid.is_empty)]

    # export as shapefile
    grid.to_file(out_vector)

    return

In order to test it I created a fake shape:
nb_cell = 6
grid_size = 30

gdf = gpd.GeoDataFrame({'geometry': [sg.Point(0,0)]}, crs='EPSG:3857').buffer(nb_cell*grid_size, cap_style=3)
dst = tmp_dir.joinpath('shapes.shp')
gdf.to_file(dst)

but when I launch my function grid(dst, 'output.shp', grid_size) I get the following error :

/usr/local/lib/python3.6/dist-packages/geopandas/base.py:39: UserWarning: The indices of the two GeoSeries are different.

Does someone see what I have done wrong?
edit
The objective is to create grid over country geometries. These geometries can include islands so multiple polygons. The idea is to create a grid on the full extends of the initial gdf and then cut all the grid cells that goes outside the country boundaries.
Based on the documentation, intersection seems to be the way to go


texts: Set background color and opacity for a WKT layer vector in OpenLayers3

texts:
I have this function:

function drawLayer(wktRegion, name,colour = "transparent"){
      var format = new ol.format.WKT();
      var feature = format.readFeature(wktRegion);
      feature.getGeometry().transform('EPSG:4326', 'EPSG:3857');

      var myStyle = new ol.style.Style({
                   fill: new ol.style.Stroke({
                      color: 'red'
                   })
                });

      var layerVector = new ol.layer.Vector({
            source: new ol.source.Vector({
            features: [feature],
            style: myStyle
        })
      });
      layerVector.set('name',name);
      map.addLayer(layerVector);
}

wktRegion can be:

POLYGON(.....)
POINT (LONG LAT)

This function give the WKT description and draw the vector layer on the map.

I am trying to set a background color for the vector layer and an opacity, but with this code I obtain always transparent vector layer.

What I am doing wrong?


texts:
I have seen over the internet a lot of examples of using Features paste them into the Vector layer and after that insert the layer to the map .

basic example:

 <script>
  var iconFeature = new ol.Feature({
    geometry: new ol.geom.Point([0, 0]),
    name: 'Null Island',
    population: 4000,
    rainfall: 500
  });

  var iconStyle = new ol.style.Style({
    image: new ol.style.Icon(/** @type {olx.style.IconOptions} */ ({
      anchor: [0.5, 46],
      anchorXUnits: 'fraction',
      anchorYUnits: 'pixels',
      opacity: 0.75,
      src: 'data/icon.png'
    }))
  });

  iconFeature.setStyle(iconStyle);

  var vectorSource = new ol.source.Vector({
    features: [iconFeature]
  });

  var vectorLayer = new ol.layer.Vector({
    source: vectorSource
  });

  var rasterLayer = new ol.layer.Tile({
    source: new ol.source.TileJSON({
      url: 'http://api.tiles.mapbox.com/v3/mapbox.geography-class.jsonp',
      crossOrigin: ''
    })
  });

  var map = new ol.Map({
    layers: [rasterLayer, vectorLayer],
    target: document.getElementById('map'),
    view: new ol.View({
      center: [0, 0],
      zoom: 3
    })
  });

In all of those examples they Feature created and put into Vector layer before we create the map but if I want in real time to add some Features.. after I created the map ? how can I do it ?


texts: How to group insert statments into transactions when using ogr PGDump driver?

texts:
I tried the -gt option to group insert statements into transactions, but apparently it doesn't make a difference:
ogr2ogr -f PGDump test1.sql test_6350.gpkg -lco SRID=6350 -lco DIM=2 -nlt PROMOTE_TO_MULTI -gt 50000
ogr2ogr -f PGDump test2.sql test_6350.gpkg -lco SRID=6350 -lco DIM=2 -nlt PROMOTE_TO_MULTI -gt 100000
diff test1.sql test2.sql is empty.
There are as many INSERT statements as the number of features in  test_6350.gpkg, ~2 million. I could grep only one line with a COMMIT in it. So next when I try psql "connection string ..." -f test1.sql how is it going to know how inserts are grouped?


texts: Same equal area projection results in different areas when shifting the center of projection

texts:
I found that using an equal area projection results in different areas for the same polygons if you shift the center of projection (ie to minimize the distortion of your area of interest when you project it).

How is it possible? This defies my understanding of an equal area projection keeping areas intact.

Example:

library(sf)
options(scipen = 20)

south_america <- st_read("https://raw.githubusercontent.com/codeforamerica/click_that_hood/master/public/data/south-america.geojson")

# Transform to Lambert azimuthal equal-area projection, default center (0° N 0° E)
south_america_laea <- st_transform(south_america, "+proj=laea +nodefs")
plot(st_geometry(south_america_laea))

st_area(south_america_laea)

  Units: [m^2]

  83375755616 8467140414780  737655036685 2782483050075 1086627258107 1134458903620  254527695525

  211024005330 1290996983204  399946127945  144935413395  177467272516
  910435563567   11429672354

# Same Lambert azimuthal equal-area projection, just shift the center of projection to center it on the continent

south_america_laea_2 <- st_transform(south_america, "+proj=laea +lat_0=-32 +lon_0=-60")

plot(st_geometry(south_america_laea_2))

st_area(south_america_laea_2)

  Units: [m^2]

  83372539973 8467111444577  737677003441 2782466797426 1086627986050 1134453525235  254528299278

  211021681369 1291025934349  399948128115  144931686414  177467118512  910429865771   11429872564

Now the areas are different! Not even a linear transformation, just slightly different numbers


texts:
How does google or other map providers with routing capabilities determine when you are not staying at the right course?
I am trying to figure out how to implement this feature in my gps that I am developing.
Right now I have the polyline with all the coordinates and the location of the device. Is there any general algorithm I can use?


texts: In WMS, is there a way to provide appropriate resolutions/zoom levels?

texts:
I am looking for a way to tell a wms client which zoom levels or resolutions that are appropriate for a layer. Something like

...
<layer>
  ...
  <resolutions>
     <resolution>1.0<resolution>
     <resolution>2.0<resolution>
     <resolution>4.0<resolution>
  </resolutions>

One "less preferred" way would be the VendorSpecificCapabilities in the WMS GetCapabilities 1.1.1 DTD, but it is dropped in the 1.3.0 XSD.

Is there another way to do this within the WMS protocol?


texts: Zero filling attributes using ArcPy/Python?

texts:
I am just learning python and I need to convert a 2 digit numeric attribute in a table to a 2 digit zero-filled attribute in a script.

For example, 9 needs to be 09.


texts: How do you use GIS to detect percentage of Land Use Change

texts:
I am a student in GIS I have a dissertation about land use change. I have questions  about using ArcMap 10 for this aim, how can using GIS to detect the percentage of land use change. For example  we want to identify land use change in Sheffield city between 1990 until 2000
Land use of Urban %30 1990 also it is changed to %40 in 2000 how can detect this result I mean (%30 and %40).


texts: Wrapping a GeoTIFF over the antimeridian in QGIS

texts:
I want to make a map of marine sampling locations for a study, and in this context show some bathymetric profiles. So in my map I have shapefiles from NaturalEarth as well as GeoTIFF bathymetry data. One location I am interested in is the Chatham Rise, and its location unfortunately extends over the edge of the map in my case:

Can I somehow achieve, that the out of bounds part of the image is "wrapped around to the other side" of the map? The left edge of the map shows my actual sampling location:

Or will I have to specify a different projection to achieve what I am trying to do? At the moment I have no projection selected at all.


texts: QGIS 3.10.0 SCP SSL Errors Downloading Landsat Mac OSX 10.14 (Mojave)

texts:
Semi-Automatic Classification Plugin (SCP) throws an SSL error when trying to download or preview Landsat images.

Looking for a workaround, this has been addressed for Windows users but not Mac.

https://semiautomaticclassificationmanual-v5.readthedocs.io/en/latest/faq.html#error-56-ssl-connection-error-unable-to-download-sentinel-2-images-why

I would imagine it is related to the Python installation that is bundled with QGIS and security settings on Mojave, it may be an issue that has been solved for other plugins.

  ERROR exception:
  url:https://ims.cr.usgs.gov/browse/landsat_8_c1/2019/019/030/LC08_L1TP_019030_20191011_20191011_01_RT.jpg


texts: Why does Contours plugin crash QGIS 2.18?

texts:
I am running QGIS 2.18 on a mac and trying to use the Contours plugin. I've never had problems with it before but then I updated QGIS using Kyngchaos.  Now Contours crashes QGIS when I choice "Filled Contours" or "Both".  It works fine for generating just contour lines.  When it crashes, it doesn't throw up an error window; the whole program just quits.

I tried stepping back to an older version of QGIS but that not only didn't solve the problem, it also created a host of other issues. I've also uninstalled and reinstalled the plugin.

How can I get Contours plugin fully functional again?


texts: Possible to use wordwrap in the identify features tool in QGIS

texts:
I was wondering if it is possible to change the appearance of the identify features tool in QGIS by applying wordwrap.

I have a field ("description") with multiple lines of information (up to 1000 characters) about archeological monuments. I would like to see these lines of text capped at 50 caracters with the reminder of the text continuing on the following lines whenever I click on an archeological monument.

Is this possible in QGIS?


texts: Creating an oblique Mercator projection using proj and QGIS

texts:
I am trying to create a coordinate system within QGIS that's aligned with a great circle between two points.
The points defining the great circle are (-44.0597,168.7192) and (-42.7365,171.4428) and, ideally, I'd like the first of these to be correspond to the origin of the new system. From what I can tell from proj.org, the following definition should work:
+proj=omerc +lat_0=-44.0597 +lonc=168.7192 +lat_1=-44.0597 +lon_1=168.7192 +lat_2=-42.7365 +lon_2=171.4428 +units=km

However, when I test this, I'm told that (-44.0597,168.7192) maps to (0,0) [correct] whereas (-42.7365,171.4428) maps to (143,223) [incorrect].
This means that the new coordinate system has the origin I wanted but is not correctly aligned with the great circle. I've also tried specifying the new coordinate system using a center point and either +alpha=-33 or +gamma=-33 (since the great circle has an azimuth from the first point of about 57 deg.) but this isn't working correctly either.
I'm sure there's something obvious I'm missing but after going round and round in circles, I'm at a loss to spot it. I'm using QGIS 3.10.2 and PROJ 6.3.0.


texts: Incorporating turn restrictions shapefile into network dataset in ArcGIS Desktop?

texts:
I've started creating routes utilizing Network Analyst with a municipal centerline. I quickly noticed that a few of the routes were incorrect as per the restrictions for turns at given intersections.

I found a municipal turn restrictions shapefile that contains what I need:
Turn Restrictions

Here is the geodatabase feature dataset which contains my centerline that I used to initially start creating routes. Elevations, speed, and one-way streets are already incorporated into this file. There is also a blank turn feature class.
geodatabase feature dataset

I have not been able to migrate the data from the shapefile into a workable format that the turn feature class can recognize.

Following the steps on this ESRI site:
Link

address what I'm dealing with concerning the shapefile turn data leads me to a dead end, notably that the second step utilizing the "turn table to turn feature class" tool requires ArcVIEW or ArcINFO table format giving me the error:

My other option was to try to load the table data from the shapefile into the turn feature class which showed that parts of the data fields match, but the shapefile (whether it is a feature class or not, is missing all the "edge" field parameters while matching with all of its remaining fields which can be seen here:

In this case, is there a step I have missed or is there a workflow within Catalogue to build this edge data from the shapefile data?


texts: QGIS expression on layer that refers to intersects on same layer

texts:
Let's say I have the scenario where

Polygons A and B intersect polygon c
There are also polygons X,Y and Z which do not intersect c.
Polygon C has an attribute value 'Priority' which I need the value propagated onto intersecting polygons (that is, onto Poly A and B).

All the data is in a single layer in a geopackage table.
I have no drama's selecting polygons A and B, I'm looking at trying to update them using the Field Calculator and an expression. But I'm having an absolute brain fade as to how to refer to self.
That is - how do I find intersecting polygons on the layer that I'm already using and get their attribute values?
My logic is to write a function and apply it to a selection set (which I would select outside of the function). So far I have the following:
def get_value_func():
# Firstly, get the set of polys I want to grab values for. Lets just use the 'c' poly for now
geometry_col = collect($geometry, , filter:="layer_name" = 'c')

#then lets test if the current geometry intersects anything within the geometry collection
if(intersect($geometry, geometry_col), update_val = <THIS IS THE BIT IM NOT SURE ABOUT>))

return update_val

Its possible I'm overcomplicating things and am happy to hear other ideas about how to achieve this. Maybe collection isn't the right data structure to be using.


texts: Trace tool being very erratic in ArcMap?

texts: Right now I am working on a project where I have a road network, and I have a layer with points on a bus route. I am trying to trace these points along the road network to create bus routes. It was working fine for the first few routes I did, but then when I started another route, it started acting really weird. It won't really go along the road network anymore, instead creating erratic lines that don't go along anything, often lines that don't even start where I originally click. The coordinate systems of all my layers are the same. What's happening here?


texts:
I have used ogr2ogr to import various formats to PostgreSQL and have used GDAL utilities. Is there a complete Java port for these tools or are they simply wrappers calling JNI to C++ code?

I have seen some Java stuff by Googling but not sure how much are they current and kept up to date.


texts: Road distance from origin to specified destination using ArcGIS network analyst?

texts:
I am looking for a way to compute network distance from a given list of points to their predefined destinations. Please note, I am NOT looking for closest distance but the distance to a destination that I have specified for each point.

Is this possible using the ArcGIS network analyst extension?


texts:
I use PostGIS, GeoServer, and OpenLayers for GIS development.

I simply reach from GeoServer to PostGIS table or view and show it with OpenLayers, but I need some filtering operation on my spatial data. For example, I have historical spatial data and filter them to start_date and finish_date like:

select .......
from ......
where date between :start_date and :finish_date

The problem is start_date and finish_date are dynamic.  User changes them from web my open layer web project.

How can I do it using GeoServer?


texts: how to include part of a map extent in layout view

texts: I have a map of Scotland and within the map layout I'd like to include Shetland within a box at the top so that I don't have to minimise the full extent of the map - I was wondering what the best way to do this would be. Previously i created a duplicate layer in the data view and zoomed it to only fit in Shetland in the layout view but this seemed a clunky way of doing this? Anyone have any advice?


texts: How much time is required to get my osm data updated

texts: I updated some roads in osm map. I updated using JOSM. But the updates are not visible in the map after 2 days. I am new to osm maps. Can anyone please tell me the possible reasons.


texts: Batch processing_Distance to nearest hub

texts: How do I go about batch processing for distance to the nearest hub? I have about 10 variables/covariates. Have tried severally and not seems to succeed.


texts: PostgresSQL/PostGIS: how can others access the data?

texts:
I work in an office of 10 people. We all use QGIS. I also use ArcGIS.

Currently all our GIS data are stored (thus duplicated) on everyone's induvidual hard drives. Any major updates I make are copied over manually from my machine to theirs.

We want to move to having the data in a single location. I thought that by installing PostgresSQL with PostGIS on my machine, and importing the data using pgAdmin, I'd be able to access the same data on my colleague's machine.
However, I can only connect to the database (in QGIS) on my machine.

My question is: how can I share my PostGIS database with others in my office?


texts: How to create geometry with geometry functions?

texts:
I see there are functions in QGIS that supposed to generate geometry, like geom_from_wkt. How do I actually display it?

When I create a new field (in a field calculator) with a formula, say

bounds(make_line(make_point($x-1,$y-1),make_point($x+1,$y+1)))

its value becomes NULL no matter what type I select.

And If I say Update existing field" and choose <geometry> in the dropdown -- nothing changes, I still do not see rectangles.


texts: What the default unit in QGIS 2.18 $area is calculated and etc?

texts:
I'm a beginner  in QGIS . We are using 2.18 version.

What are the basic units for area and slope, etc? SI, so means m, m^2, °? Or I can set it by myself? If it's possible to set it, where?


texts: Optimizing a nearest neighbor search in PostGIS for a large dataset

texts:
The below code is working well for smaller datasets, however I want to use this with a 30 GB table. Every time I run it using the larger table I run out of SSD space, the table gets to 127 GB in size!

Basically what I am try to achieve is calculate the minimum distance from each land parcel and the closest bus stop. I then want to save that distance in a separate column. Land parcels that do not have a bus stop within 2km, should have its distance set to -9999.

CREATE TABLE parcels_all_shapefile_bs AS
SELECT a.ogc_fid,
       a.poly_id,
       a.title_no,
       a.insert,
       a.update,
       a.rec_status,
       a.wkb_geometry,
       MIN(COALESCE(ST_Distance(a.wkb_geometry, b.wkb_geometry), -9999)) AS bs
FROM parcels_all_shapefile_cp a
    LEFT OUTER JOIN filter_busstops b ON
        ST_DWithin(a.wkb_geometry, b.wkb_geometry, 2000)
GROUP BY a.ogc_fid,
       a.poly_id,
       a.title_no,
       a.insert,
       a.update,
       a.rec_status,
       a.wkb_geometry;

Can anyone see what I am doing wrong and how I could optimize the above code?


texts: Changing Layer CRS Causes Misalignment and Scaling Issues

texts:
Update:
My final solution was to re-download the data sets and all my alignment issues were solved.
Background:
I am working on a personal project overlaying the location of churches/parishes on various counties in New Jersey. The parishes are contained in a .geojson file and the state counties come from a shapefile provided by the state.
Problem:
When both layers are loaded the parish layer is way up to the left of the counties and both have different CRSs. The parish is EPSG:4326 - WGS 84 while the counties are the same as the project crs EPSG:3857 - WGS 84 / Pseudo Mercator. When I change the crs of the parish layer to the projects there are alignment issues and scale issues. The parishes are much closer (but still not over the counties) and to the bottom left. Further the scaling goes wrong and all the data points are much to close to each other.
How do I resolve this alignment issue?
The final result is I will be making a map to print out, which I can do now my manual exporting the layers at a fixed scale. Then I can manually align the two exports in my image editing software. I would prefer to figure out this issue mainly to learn more about these things.
Key Information and Links:
Parish Layer CRS: EPSG:4326 - WGS 84
Project CRS: EPSG:3857 - WGS 84 / Pseudo Mercator
Base data source comes from State of New Jersey - GIS Open Data State Counties.
Link Here:
https://njogis-newjersey.opendata.arcgis.com/datasets/5f45e1ece6e14ef5866974a7b57d3b95_1/geoservice?geometry=-78.816%2C39.583%2C-70.472%2C41.049
I made the parish data source using geojson.io. You can view the data source one either:
GitHub: https://github.com/pulsence/diocese-of-metuchen-parish-map or
GeoJson: http://geojson.io/#id=github:pulsence/diocese-of-metuchen-parish-map/blob/master/parish%20locations.geojson&map=10/40.6289/-74.7065
I am using for this project QGIS 3.0.3-Girona.


texts: I have a 6 mile line that I need to make a point every 500 ft on it.  Any idea what tool I can use to do this?  I have arcmap 10.1 standard.


texts: Assign Icons to specific points dynamically in OpenLayers

texts:
I have a json file describing features on a map (routes, points etc)

points include places of interest such as airport, harbor etc.

Is there a way to assign icons to these points?

example: have a path to an icon in my json file as an attribute to that point and request that that icon be used as a marker for that point?

I used the advice below by @Aragon, still can't get it to work..

my json looks like this:

{
"type": "FeatureCollection",
"features": [
    {
        "type": "Feature",
        "id": "id0",
        "properties": {
            "graphic": "../../icons/pushpin.png",
            "title": "blah",
            "transportation": "Blah",
            "description": "blah"
        },
        "geometry": {
            "type": "Point",
            "coordinates": [
                24.067939338686,
                35.525197669582
            ]
        },
        "crs": {
            "type": "name",
            "properties": {
                "name": "EPSG: 900913"
            }
        }
    }
]}

and my javascript looks like this:

 var context = {
        getGraphic: function(feature) {
            return feature.attributes.graphic;
        }
    };
    var template = {
        'externalGraphic': '${getGraphic}',
        'graphicWidth' :15,
        'graphicHeight':15,
        'fillColor':'#669933',
        'fillOpacity':.8,
        'strokeColor':'red',
        'strokeWidth':5,
        'pointRadius':6

    };
            var data_vector_style = new OpenLayers.Style(template, {context: context});

            // Style Map tells vector Layer what style objects to use when features are
            //  in a certain state (states can be default, select or temporary
            var data_style_map = new OpenLayers.StyleMap({'default':data_vector_style});
            var dataLayer =  new OpenLayers.Layer.Vector("My GeoJson data from server", {
                projection:  new OpenLayers.Projection("EPSG:4326"),
                styleMap : data_style_map
            });

            var features = getFeatures("../../controller?action=GET_USER_LOCATIONS");
            //alert(features);
            dataLayer.addFeatures(features);

if I change any of the other styling in template, it works fine. if I hardcode the marker in the styling    'externalGraphic': '../../icons/pushpin.png', it works fine. when I try to use attribute replacement the marker isn't used... Am I doing it wrong in the json object?


texts: JTS Geometry "Contains" not detecting point of an Inner polygon on edge of Outer Polygon

texts:
I was wondering if a JTS expert would be able to explain how to detect the below intersection of an inner polygon inside an outer polygon.
I would have thought that at least the "contains" check be false due to the point at 1000/500 being along the right-edge of Polygon1.
Polygon1: POLYGON ((0 0, 1000 0, 1000 1000, 0 1000, 0 0))

Polygon2: POLYGON ((500 500, 1000 500, 600 600, 500 600, 500 500))

Results in:
polygon1.contains(polygon2) = true

polygon1.covers(polygon2) = true

polygon1.touches(polygon2) = false

The only way I can detect this this inner polygon touching the edge of polygon 1 is by iterating through polygon2 points and creating separate points and checking:
polygon1.touches(eachPolygon2Point) but I'm worried about the performance impacts of doing this on a routine basis. Note sometimes polygon2 might be outside polygon1 so I don't really want to rely on say creating polygon2 as an inner ring of polygon1 every time either.


texts:
I have multiple S57 Ecdis files and would like to extract some data out of them. In specific some 'depth' data which is bound to lat/lon coordinates.
I will need to bulk process all the data and group them based of their locations. Basically we have a map where polygons are drawn on and for each polygon we need to calculate the min/max/median of all the depth-values that are inside it.

Do you have any recommendations on how to do this?
Would it be possible to convert it to NetCDF file or to some kind of raster?

The program will be written in Java so preferably I’d like to convert it using some Java libs.


texts: Using sym.renderer.fields with ArcPy and ArcGIS Pro gives NameError: attribute 'fields' not supported on instance of UniqueValueRenderer

texts:
I have a script that updates the symbology of layers in ArcGIS Pro 2.5.0 with Python 3.6:

if lyr.name == "Oiling Zone":
    sym = lyr.symbology
    sym.updateRenderer('UniqueValueRenderer')
    sym.renderer.fields = ['ZSURFDOO']
    for grp in sym.renderer.groups:
        for item in grp.items:
            value = item.values[0][0]
            if value == "Heavy":
                item.symbol.color = {'RGB': [250, 52, 17, 100]}
                item.symbol.size = 1.5
                item.label = str(value)
            if value == "Moderate":
                item.symbol.color = {'RGB': [255, 170, 0, 100]}
                item.symbol.size = 1.5
                item.label = str(value)
            if value == "Light":
                item.symbol.color = {'RGB': [255, 255, 0, 100]}
                item.symbol.size = 1.5
                item.label = str(value)

    lyr.symbology = sym

This code previously worked, but has started throwing this error on the 4th line:

NameError: The attribute 'fields' is not supported on this instance of UniqueValueRenderer.

Why am I getting this NameError?

Documentation says fields is an attribute of UniqueValueRenderer.  This is how the layer is displayed in the TOC:


texts: commuteDistance in gdistance package in R giving error

texts:
I am currently trying to run a function that "calculates the expected random-walk commute time between nodes in a graph. It is defined as the effective distance (resistance distance) between the selected nodes multiplied by the volume of the graph, which is the sum of the conductance weights of all the edges in the graph (Chandra et al. 1997). The result represents the average number of steps that is needed to commute between the nodes during a random walk"
To run commuteDistance from the gDistance package, I need to create a transition object.
Step 1: I loaded my rasterlayer
> jul <- EVIrec[[6]]
> jul
class       : RasterLayer
dimensions  : 642, 382, 245244  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : 461951, 843951, 892583.3, 1534583  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Jul
values      : 1, 3  (min, max)

Step 2: I created a 'transition object' as mentioned here
> tr <- transition(jul, mean, directions = 8)
> tr
 class       : TransitionLayer
 dimensions  : 642, 382, 245244  (nrow, ncol, ncell)
 resolution  : 1000, 1000  (x, y)
 extent      : 461951, 843951, 892583.3, 1534583  (xmin, xmax, ymin, ymax)
 coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
 values      : conductance
 matrix class: dsCMatrix

Step 3: I used a geographic correction approach on the transition layer. Details here
> julcorr <- geoCorrection(tr, type = "r", multpl=F)
> julcorr
 class       : TransitionLayer
 dimensions  : 642, 382, 245244  (nrow, ncol, ncell)
 resolution  : 1000, 1000  (x, y)
 extent      : 461951, 843951, 892583.3, 1534583  (xmin, xmax, ymin, ymax)
 coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
values      : conductance
matrix class: dsCMatrix

I then used the commuteDistance function as mentioned here and I get this error as shown below.
> jul_comm <- commuteDistance(julcorr, centroids_jul)
Error in LU.dgC(a) : cs_lu(A) failed: near-singular A (or out of memory)
In addition: Warning message:
In .rD(x, coords) :
  46 out of 47 locations were found in the fully connected transition matrix. NAs introduced.

Ultimately, I would like to run a function that does an operation similar to CIRCUITSCAPE, but within R and then use these results for further analyses.


texts:
I have created a model that (I want) to take a point, selects the polygon it is within, selects a different point that is within the same polygon that contains the data I am after.  I then want to join the tables and get rid of unnecessary columns.  I have toggled the Field Info checkboxes in the Make Feature Layer tool, however, every time I run the model, those boxes reset.  I am assuming this is because the tables that I am trying to join and make into a feature layer are temporary, so when the Make Feature Tool is ran it is going to its default settings as it is looking at a new table every time.

How do I make the tables I am ultimately joining permanent?


texts: Get km route between two points with lat / lng

texts:
I'm doing the following query to get the distance between two points with pgrouting:

SELECT
        sum(cost) AS cost
FROM pgr_dijkstra(
                'SELECT id, source, target, km as cost FROM at_2po_4pgr',
                50347, 32988, false, false
        );

It turns out that the source and destination comes from nominatim geocode coordinates in, and I do not know how to use the coordinates as parameter in sql.

Can anyone help me?


texts:
The WKT specification allows both period and comma to be the decimal separator in floating-point numbers (section 7.2.1):

  <decimal point> ::=  <period> | <comma>

However, a comma could result in ambiguous WKT representations.
In this example, the polygon is specified with five points, but the WKT could also be parsed as having three points:

POLYGON((0 0,0 1,1 1,1 0,0 0))

I'm assuming that nobody deliberately uses , as decimal separator, but is that really true?

Is there any reason to try to parse numbers in WKT with ,?
Or shouly I assume that the specification is just buggy?


texts: Displaying/feeding US census/tiger data into Google Maps API?

texts:
Let's say I want to display age, income, etc based on census block groups on Google maps.

How can I do this and also to allow simple queries?


texts: Finding optimal number of clusters for wekaKmeans unsupervised classification in the Google Earth Engine

texts:
I want to find the optimal number of clusters for wekaKmeans unsupervised classification in the Google Earth Engine. The reason is that I want to see which cluster division is the best.
Here is my code
// Create a collection by giving a list to the constructor.
var image_comb = ee.ImageCollection([image, image2,image3,image4,image6,image7,image8,image9,image10,image11,image12,image13,image14,image15]);

print(image_comb);

// Compute a mean image and display.
var fuzzyOverlay = image_comb.mean();

print(fuzzyOverlay);

Map.setCenter(-88.24153442382814, 40.11773641362424, 10);
Map.addLayer(fuzzyOverlay.clip(champaign_border), {min:0,max:1, palette:['blue','green','red']}, 'fuzzyOverlay');
// Create Unsupervised classifaction into mean  image

var input = fuzzyOverlay;

// Make the training dataset.
var training = input.sample({
  scale: 10,
  region: champaign_border,
  numPixels:50000
});

// Instantiate the clusterer and train it.
var clusterer = ee.Clusterer.wekaKMeans(3).train(training);

// Cluster the input using the trained clusterer.
var result = input.cluster(clusterer);
var result = result.clip(champaign_border);

// Smooth with a mode filter.
var postcl = result.focal_mode();
// Use erosion and dilation to get rid of small area.
var image5 = postcl.focal_max(4).focal_min(6).focal_max(4);
var image5 = image5 .clip(champaign_border);
// Export the image, specifying scale and region.
  Export.image.toDrive({
  image: image5,
  description: 'Unsupervised_fuzzyOverlay',
   scale: 10,

  region: champaign_border
  });


texts: How to label raster cells using ArcGIS Desktop or QGIS?

texts:
I want to generate a small sample raster (eg. euclidean distance or cost distance) to explain the method of cost-distance analysis, and I'd like to label the pixels with the corresponding cost-value.

I already tried to lay a fishnet over the raster in ArcGIS, but I didn't find a good way to copy the values into the attribute table.


texts:
I'm using Polygon to Raster tool to create a raster (10x10). The problem is that some cells of the output raster is not assigned any value.

My understanding is that the tool always uses the cell center to decide the value of a raster pixel.
Attached is the result when I overlay both the output raster and Input polygon on ArcMap. (Black is the NoData cell).

Anyway to include all those cells that fall within the polygon when generating the raster?
I tried to use Feature Vertices to Points to generate the point and convert the points to raster, but the vertex of the polygon that I used only have 1 vertex.

Another item to note is when using smaller cell size (1x1), the raster cell generated are good but when I resample it to 10x10, the output raster from resampling is the same as above raster.


texts: Large Scale Geocoding and Processing in ESRI

texts:
Ok, so I guess this kind of a informal query/survey about how big a datasets you are using in your ESRI worlds...
I am building and maintaining a statewide dataset, where I have to process down to the individual house level, not parcel level but multiple mailing addresses per parcel for our systems. In many places I am using theoretical addresses calculated from street network or USPS AMS/AIS data. So my Address List is roughly 13.5 million addresses and growning monthly or quarterly.
Is anyone out there right now maintaining a live system of address/properly lookup information that is this large in a continuous dataset? I would love to collaborate or talk more about how others are handling such a large dataset. I am seeing issues where ESRI software seems to be blowing up when I try to perform tasks such as intersects or spatial joins. ESRI says they don't see these kinds of issues but I have had these issues since back to 9.3.1 so I can't be the first/only person doing this since I can recreate it across multiple machines. My Platform right now is ESRI ArcGIS 10 on the Desktop, talking to ArcSDE 9.3.1-sp1 on a SQL2008 backend using the GEOMETRY spatial object. So I am not doing anything really exotic; but still seems to me that in some areas I maybe am pushing the envelope.[Further] What I am interested in know is what are other people doing to optimize there processes for dealing with these datasets. I am going to be adding upwords of a million records a month going forward, and while Geocoding etc isn't a problem when you start running other processes and linking data for further analysis you start dealing with complex joins. Well, you output data from Intersects/Overlays/Identities using Only_FID and you get a thin middle table to join too; but when you start trying to divide and conquer the creation of that table you start to hit issues where you need to divide your source data into working areas but then you have repeating IDS that you can't merge back; so you are left with smaller blocks of data that you can't easily make whole again. Thinking about options that break the data down to County-by-County scale, then using spatial views to join it back together etc... Just curious if other users are looking at the same kinds of problems on such a large scale but on small footprints.


texts: CartoDB merge tables (column join) returns import error

texts:
I have been trying to merge two tables using the column join in Carto DB.

The first table is a table of US counties that contains the polygon geometry and the FIPS id code. The second table contains the FIPS id as well as information about population density. The common column "FIPS" is string in both tables. There are just over 3,000 rows in each table.

When I try to merge them i CartoDB it appears to begin executing successfully, then it takes me to a page that says "What a pity, the page you are looking for doesn’t exist." I am able to click back to my dashboard, and a status bar still shows that "creating table" is executing.

After some time a final error appears as: "Oops! An error occured while importing your file."

Has anyone experienced this? Any ideas on what is causing the error and how to get around it?


texts: How to find max value for id appearing anywhere in another column?

texts:
I hope that title is clear. If not, please change.

I have table, which look like this:

I would like to create a select query, where if gid is in anywhere in ancestor column, select max level from this rows.
For example:
Gid is 727 is four times in ancestor column. Max level from this rows is 3, so the value in the new column, let's says level_a, where gid is 727, will be 3.
For gid 725 the same but for gid 717 will be 0.

I tried to use unnest, join, subquery but nothing worked.
I have read about lateral join but I don't have any idea, how I can use it.

I'm using PostgreSQL 10 with PostGIS 2.5.2


texts: Importing DXF - Polylines with elevation data

texts:
When importing a DXF file (containing 3D contour data) to QGIS, none of the various methods described elsewhere are retaining elevation data as an attribute.

e.g. DXF2SHP, GRASS v.in.dxf, opening natively as a vector file in QGIS etc. DXF2XYZ is not applicable as the file does not contain point data.

The DXF file contains polylines with elevation data only. All line objects are imported successfuly and there are no errors.  Are there any other utilities?


texts: How to get value and store it in a variable from a line edit or text edit in pyQT?

texts: How to get value and store it in a variable (in main python plugin) from a line edit or text edit box that were created in QT designer ?


texts: Google Earth Engine, get a list of unique values in a list

texts:
I'm trying to figure out a way to get a list of unique values in a ee.List() object. For example:

input=ee.List([1,1,1,4,4,5,6,6])
output=[1,4,5,6]

I originally thought I could use the ee.List.distinct() function, but this does not appear to be an actual valid function. Since then I've tried a variety of different approaches including a variation of the solution to this related question to no avail.

My current approach involves an iterator function that appends newly encountered values to a new "unique values" list, however it just seems to be recognizing every value as unique and outputting the original list.

Here is my code:

var test_list=ee.List([1,1,1,4,5,6,6])
var iterlist=ee.List([])

function unique_list(inputItem, distinct_values){
  distinct_values=ee.List(distinct_values);
  var number=ee.Number(inputItem);
  if(distinct_values.contains(number)===true){
    distinct_values=distinct_values
  }
  else{
    distinct_values=distinct_values.add(number)
  }
  return distinct_values;
}

var result=test_list.iterate(unique_list, iterlist)
print('unique list is', result)

Any help finding a solution to this seemingly simple problem would be greatly appreciated.


texts: Citing background layer from ArcGIS Online that is used in ArcMap

texts: How do you cite a background map layer that is used in ArcMap that is retrieved from ArcGIS Online?  For example the background layer is suppose to be credited from the U.S. Geological Survey and in the description it was digitized from the U.S. Geological Survey Professional Paper 1183.


texts:
I try to using pgrouting to find route, first I'm using assign_vertex_id and when I calculate the shortest route, i didn't find the route in my pgAdmin data output, and in the messages :
"row number -1 is out of range 0..-1
Total query runtime: 15 ms.
0 rows retrieved."

the code:
SELECT * FROM shortest_path('
   SELECT gid AS id,
          source::int4 AS source,
          target::int4 AS target,
          length::float8 AS cost
   FROM jalan', 1,10, false,false);

and then i didn't use assign_vertex_id but there's an error,
(i'm following step on), the code:
CREATE TABLE node AS
   SELECT row_number()OVER (ORDER BY foo.p)::integer AS id,
          foo.p AS the_geom
   FROM (
      SELECT DISTINCT road_ext.startpoint AS p FROM road_ext
      UNION
      SELECT DISTINCT road_ext.endpoint AS p FROM road_ext
   ) foo
   GROUP BY foo.p;

error: "ERROR: syntax error at or near "OVER"


texts: The simplest way to uninstall any and all ArcGIS products?

texts:
What it is the simplest and most straightforward way to uninstall any and all ArcGIS products more or less in one go that might be present on a Windows machine? perhaps scriptable?

Add/Remove Programs from the Control Panel is straightforward, but when there a several products installed it's a process of click-click-click, wait for product 1 to be uninstalled, which could take some time, then click-click-click for #2, and so on. This is okay for one machine but when there are several it's real pain.


texts: QGIS: Importing DXF/DWG file without losing extra attribute information

texts:
Let me be clear here I know how to import a DWG/DXF file here the problem lies that theren is extra data that gets lost. One layer in petuclar is what we want to import the Layer XXX and separate into four different shapefiles based on the EVAL parameter. The problem is that when we import the dxf file with the features, the extra data gets lost, it is not available in the attribute table.

Picture of the data:

I have all object selected here to "Censor" the data.

How do we solve this?


texts: Optimizing Field Update - Update Cursor, ArcPy

texts:
I am trying to update a field in a (point) feature class (e.g. subLayer) using attributes and spatial proximity from another (line) feature class. For this, I'm currently leveraging arcpy tools and an update cursor. This is currently taking over 30 minutes to run. The subLayer feature class contains 7000 records and transLayer contains 7500 records; so these aren't massive feature classes, though they are continually growing as they are developed. Any suggestions on how to optimize this code?

import arcpy
from arcpy import env
import time
start = time.time()

arcpy.env.OverwriteOutput = True

defaultGdbPath = 'C:\Topo_Check_V5.gdb'

subLayer='C:\Topo_Check_V5.gdb\Subs' #point feature class
transLayer='C:\Topo_Check_V5.gdb\TLS' # line feature class

#add in Line_Count field
arcpy.AddField_management(subLayer, "Line_Count", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

#TLineCountField = "Line_Count"

#where clause
where_13_8 = ' "VT" = 13.8 '

#MakeFeatureLayer

arcpy.MakeFeatureLayer_management(transLayer, "transLayer_lyr")
arcpy.SelectLayerByLocation_management("transLayer_lyr", 'WITHIN_A_DISTANCE', "subLayer_lyr", .000002, "NEW_SELECTION")
arcpy.SelectLayerByAttribute_management ("transLayer_lyr", "SUBSET_SELECTION", where_13_8)

with arcpy.da.UpdateCursor (subLayer, ["Line_Count"]) as LineCountList:
    for subrow in LineCountList:
        result=int(arcpy.GetCount_management("transLayer_lyr").getOutput(0))
        #print result
        subrow[0] = result
        LineCountList.updateRow(subrow)
print "Line_Count Updated"
del LineCountList

end = time.time()
print (end - start)/60

In an effort to speed this up, I got a little ahead of myself and took advice to bring the selection statements out of the for loop. Although this did speed up the processing time, I ended up losing the functionality of the original script. I do want to update the unique count of all selected transLayer features (here based on proximity, but for other blocks also based on proximity and attribute data) unique to each record in the the sublayer feature class.

Here is the script that I would like to optimize. Maybe using a spatial join route rather than an update cursor will be the most efficient way to go. I am open to deviating from arcpy if there may be a better Python alternative. I am also curious if dictionaries could help speed things up.

import arcpy
from arcpy import env
import time
start = time.time()

arcpy.env.OverwriteOutput = True

defaultGdbPath = 'C:\Topo_Check_V5.gdb'

subLayer='C:\Topo_Check_V5.gdb\Subs' #point feature class
transLayer='C:\Topo_Check_V5.gdb\TLS' # line feature class
ppLayer='C:\Topo_Check_V5.gdb\PPS' #point feature class

###### Line Counts: Count of connecting lines within distance of point feature class ################

#add in Line_Count field
arcpy.AddField_management(subLayer, "Line_Count", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

TLineCountField = "Line_Count"

#MakeFeatureLayer
arcpy.MakeFeatureLayer_management(transLayer, "transLayer_lyr")

with arcpy.da.UpdateCursor (subLayer, [TLineCountField, "SHAPE@"]) as LineCountList:
    for subrow in LineCountList:
        arcpy.SelectLayerByLocation_management("transLayer_lyr", 'WITHIN_A_DISTANCE', subrow[1], .000002, "NEW_SELECTION")
        result=int(arcpy.GetCount_management("transLayer_lyr").getOutput(0))
        #print result
        subrow[0] = result
        LineCountList.updateRow(subrow)
print "Line_Count Updated"
del LineCountList


texts: How to re-reproject coordinates that are coming from GeoServer

texts:
I have a GeoJSON as follows:

  "type":"FeatureCollection","totalFeatures":131445,"features":[{"type":"Feature","id":"2009BuildingsCOS.1","geometry":{"type":"MultiPolygon","coordinates":[[[[1308926.9929020153,171186.6402624961],[1308925.12256642,171176.955270945],[1308923.6956225059,171177.2343108883],[1308918.0723277144,171148.12967222783],[1308875.0080097103,171156.45018241392],[1308880.5196885085,171184.97062908477],[1308905.5063113617,171180.14093414505],[1308907.491334925,171190.41011768577],[1308926.9929020153,171186.6402624961]]]]}

Any suggestions about what coordinates I have?

Furthermore, I need to convert the coordinates to WGS (also known asEPSG 4326`). Ideally, I want to find a way to output the coordinates as GeoJSON directly from my GeoServer. The current GeoJSON with the unknown coordinates is here:

http://ec2-54-69-8-151.us-west-2.compute.amazonaws.com:8080/geoserver/WRIA9/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=WRIA9:2009BuildingsCOS&maxFeatures=50&outputFormat=application/json


texts: Is it possible to Automate label settings with ArcMap 9.3.1

texts:
We have several users that are using ArcMap, and have asked me to create a fairly complicated label scheme.  I have created several classes, each with different Expressions describing the labels.

I would like to automate the steps of creating all of the classes, and setting all of the Expressions.

Is there a way to do this with ArcMap 9.3.1.

----Updated----

Right now, to deploy, I am going to send each user a bunch of files, with some instructions.  The instructions will look like:

Instructions to change labeling for Layer1

From ArcMap TOC, click "Display" tab
Open Layer Properties for Layer1
Go to the "Labels" tab
Change the "Method" dropdown to "Define classes of features and label each class differently"
Click "Add..." to add a new class
Enter "ClassA" as a name, and press OK
Click "Sql Query..." to define the class
Click "Load..." file and select the "ClassA.exp" file
Press OK, to close the SQL Query dialog
Click "Expression..." to define what is displayed as a label
Click "Load..." and select the "ClassA.lxp" file
Click OK to close the Label Expression dialog

Repeat these steps for ClassB, ClassC, etc...

Is there anyway to make this easier for the users?


texts:
We have several DWG files that we would like to use as Layers in GeoServer.

Therefore, what is needed is to specify the coordinate of a certain point, rotate and scale the drawing prior to exporting it as SHP.

By clicking on Tools and Geographic Location, a certain point on the drawing can be geo-referenced.

Now, how can the other two tasks best be achieved?

Thank you!


texts: How to work with GeoJSON, Leaflet and pixel coordinates?

texts:
I'm a noob trying to piece together a - as I thought - simple thing:
A Leaflet map based on a raster picture with markers coming from a GeoJSON file. As this is a fictional map, I use the CRS.Simple in Leaflet and thus just use my picture pixels as coordinate system.
From some examples I put together the following:

    var map = L.map('map', {
        crs: L.CRS.Simple,
        minZoom: -1
    });

//- set the image width and height to w and h
    var w = 1700;
    var h = 2200;

    //the stuff to be able to work with xy instead of latlon
    var yx = L.latLng;
    var xy = function(x, y) {
        if (L.Util.isArray(x)) {    // When doing xy([x, y]);
            return yx(x[1], x[0]);
        }
        return yx(y, x);  // When doing xy(x, y);
    };

    var bounds = [xy(0, 0), xy(w, h)];

    //this makes the map fit the bounds we just defined
    map.fitBounds(bounds).setMaxBounds(bounds);

// insert image name here
    var image = L.imageOverlay('img/background.png', bounds).addTo(map);

    //create map icon class
    var TaladasIcon = L.Icon.extend({
    options: {
        iconSize:     [50, 50],
        iconAnchor:   [25, 25],
        popupAnchor:  [0, -25]
    }
});

    //define icons
    var metropolisIcon = new TaladasIcon({iconUrl: 'img/metropolis.png'}),
        cityIcon = new TaladasIcon({iconUrl: 'img/city.png'}),
        townIcon = new TaladasIcon({iconUrl: 'img/town.png'}),
        villageIcon = new TaladasIcon({iconUrl: 'img/village.png'}),
        innIcon = new TaladasIcon({iconUrl: 'img/tankard.png'});

//- load geoJson file
    L.geoJson(json-test, {
        onEachFeature: layer.setIcon(metropolisIcon)
    }).addTo(map);

My GeoJSON file looks like this so far:

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {
        "name": "Cityname",
        "description": "This is a very important city",
      },
      "geometry": {
        "type": "Point",
        "coordinates": [
          719,
          1210
        ]
      },
    },
  ]
}

There is obviously something wrong, as I only see the background image but no markers from the GeoJSON. I think it might have to do with part that exchanges latlng for xy, but honestly I'm not sure what exactly happens there.
It does work if I define a marker directly in the script, but not from the GeojJSON. Which is a bummer because I plan on having quite a few markers and being able to write them down in a table on geojson.io would save me a lot of time.

Can anyone point out what I'm doing wrong here?

EDIT: Second try with AJAX, still not working:

var map = L.map('map', {
  crs: L.CRS.Simple,
  minZoom: -1
});

var w = 1700;
var h = 2200;

var yx = L.latLng;
var xy = function(x, y) {
  if (L.Util.isArray(x)) {    // When doing xy([x, y]);
    return yx(x[1], x[0]);
  }
  return yx(y, x);  // When doing xy(x, y);
};

var bounds = [xy(0, 0), xy(w, h)];
map.fitBounds(bounds).setMaxBounds(bounds);
var image = L.imageOverlay('img/kisandra.png', bounds).addTo(map);

var TaladasIcon = L.Icon.extend({
  options: {
      iconSize:     [50, 50],
      iconAnchor:   [25, 25],
      popupAnchor:  [0, -25]
  }
});

var metropolisIcon = new TaladasIcon({iconUrl: 'img/metropolis.png'}),
    cityIcon = new TaladasIcon({iconUrl: 'img/city.png'}),
    townIcon = new TaladasIcon({iconUrl: 'img/town.png'}),
    villageIcon = new TaladasIcon({iconUrl: 'img/village.png'}),
    innIcon = new TaladasIcon({iconUrl: 'img/tankard.png'});

var jsonMetropolis = new L.geoJSON();

$.ajax({
dataType: "json",
url: "json-test.geojson",
success: function(data) {
  $(data.features).each(function(key, data) {
      jsonMetropolis.addData(data);
  });
}
});

jsonMetropolis.addTo(map);


texts: I wish to generate a road with approximately 1090 XY coordinates, my only problem is that my XY to line tool will not execute, i get this message 'Error 999999: error executing function. Failure to execute (XY To Line).' what could the problem be?


texts: Finding Segment information on a Route between 2 places

texts:
Is there a way to find Segment information on a Route between 2 places using Google Maps API?

I have 2 places/addresses/latlong. I am looking for the direction between the two. For this direction, I need type information about the segment:

Is it a road
Is it a highway
Is it a bridge
...


texts: How to Import the Multipolygon with GeoJSON code in Google Earth Engine?


texts: Why does IDatasetContainer2.CreateDataset fail with ArcObjects at ArcGIS 10.1?

texts:
I have a small application used to check and generate network, created with the help of the following tutorial:
http://help.arcgis.com/en/sdk/10.0/a...00000w7000000/

It worked just fine with ArcGIS 10, but since I had to upgrade to 10.1, the application fails on this line:

INetworkDataset networkDataset = (INetworkDataset)datasetContainer2.CreateDataset(deDataset);

The error code is "System.Runtime.InteropServices.COMException was unhandled Message=Exception from HRESULT:
0x80042262", which really isn't helpful. Is there something new I could have missed, or some help on what could cause such an error ?


texts: QGIS: digitized object disappear when filling in attributes

texts:
When digitizing, for example a point, the attribute window opens for you to fill in the pre-defined columns/information. But no point is shown in the map while you do this in QGIS 2.6. It was shown in 2.2 and I wonder if there is a way to have the object shown while the object attribute dialog is open? (same goes for lines and polygons)

Picture below describes the situation where inside the red circle I have added a shape point feature and now prompted to fill out my attributes. The red circle was drawn on the picture only for illustration.

I run a network installation of QGIS 2.6.1, Win7, 32bit.


texts: How do I separate the transform from the data for gdaltransform?

texts:
I have a GeoPDF that is understood by gdaltransform. I would like to separate the transform from the data. Is it possible to convert the file to its transform only, in a format that gdaltransform can understand?
I would like to convert the PDF to SVG, and then translate screen coordinates from the SVG to latitude, longitude based on the original GeoPDF metadata.


texts: Problem with overlaying polygons onto DEM with qgis2threejs plugin

texts: I am trying to export a scene using qgis2threejs, adding a polygon layer overlay on top of a DEM and a satellite image. After exporting, it looks like the polygons don't overlay onto the DEM correctly. I get gaps like in the photo below. I have tried messing with the Z coordinates of the layer (I have it set so it's relative to the DEM). The vertical exaggeration of the DEM has no effect. Does anyone know a solution to this?


texts: How to load my localhost (MapServer) on my website?

texts:
My localhost is http://localhost:8080/
and I am using MS4W 3.0.3 on this localhost for making the application as Hello world, single map, more than one layer map, and so on.

Now my folder is ms4w in C drive like:- c/ms4w  and I want to upload this to my own server on internet, I thought it would be a simple application, and not running application as like my localhost, but in my server not have MapServer, it should be if yes how to run MapServer on my website and I want to run application as that I create in my local host.

In my local host it's working fine, because its run on MapServer on my localhost.

Hence everybody known this problem please please help me how can I solve the problem, its urgent.


texts: How to implement 1:n joins with string attribute concatenation?

texts:
I have a question concerning pyQGIS. I'm a QGIS long-time (self-learned) user, but just a pyQGIS newbie. I want to implement 1:n spatial join with attribute concatenation in the table. I have a vector points layer with multiple points on the corresponding polygons, I already made a spatial join of the points with the polygons, so I have now a point layer with the key field of the polygon (column 10 in the table).

Now I just want to search and compare the table attributes in column 10 and if they are the same, concatenate the text in column 5. Should be easy, but for me it isn't.
at first I just wanted to compare the column 10 and if it is the same - print how many time it repeats, so my first approach:

#!/usr/bin/env python
# -*- coding: cp1250 -*-

layer = qgis.utils.iface.activeLayer()
features = layer.getFeatures()
count=1
for feature in features:
    # fetch attributes
    attrs = feature.attributes()
    attrs1 = features.next()

    # attrs is a list. It contains all the attribute values of this feature
    print attrs[10]
    print attrs1[10]
    print "FID:"+str(feature.id())
    if attrs[10] == attrs1[10]:
        count+=1
        print count
    else:
        count=1
        print count

the program works, but it skips each 1 row.

the next approach (with help of Krzysztof Maślanka):

#!/usr/bin/env python
# -*- coding: cp1250 -*-
layer = qgis.utils.iface.activeLayer()
features2 = layer.getFeatures()
features = layer.getFeatures()
powt=1
for feature in features:
    attrs = feature.attributes()
    for feature2 in features2:
        # fetch attributes
        attrs1 = feature2.attributes()
        # attrs is a list. It contains all the attribute values of this feature
        if attrs[10] == attrs1[10]:
            powt+=1
            print "FID:"+str(feature.id())
            print attrs[10]
            print attrs1[10]
            print powt
        else:
            print "FID:"+str(feature.id())
            print attrs[10]
            print attrs1[10]
            powt=1
            print powt

The problem with the above program is that it should compare every row to each other, but the primary loop exits after the sub-loop is finished, so it just compares the 1st row to the table, then exits.

Can anybody help with this? It would cover the lack of this in QGIS. OR maybe somebody knows how to make a multiple 1:n spatial join with attribute concatenation in GRASS GIS or some other OS4GEO software (I didn't find it, but why should one ram against open door if it exists :-)?


texts: ArcGIS Dissolve gives Error 999999 on one computer but works on a different one

texts:
I'm trying to run the geoprocessing tool "Dissolve" on a single feature class in a File GDB, yet I'm getting the error 999999 "the operation was attempted on an empty geometry", but there are no empty geometries (they all have area greater then zero), and for some reason the same tool on the same GDB seems to work on a different computer.

I've checked any possible differences in the ArcCatalog settings or in the Environment variables and couldn't find any. Also tried repairing the feature class's geometries, splitting multipolygons into single polygons, deleting relativly small features or features with a lot of vertices. Did not help.

Anyone know what may cause this problem?


texts: Hide/Show leaflet markers based on specific values

texts:
I have a Python list of dictionaries with a number of IPs and "Scores" for each IP.

I'm adding them to my front end map / leaflet using jinja2 like so:

var map = L.map('map')

    .addLayer(mapboxTiles)
    .setView([42.3610, -71.0587], 4);

    var data = {
      "type": "FeatureCollection",
      "features": [
        {% for obj in data %}
        {
          "type": "Feature",
          "properties": {
              {% for obj2 in obj.interfaces %}
                {% for obj3 in obj2.ip_addresses %}
            "full_name": "ID: {{ obj3.value }}",
                {% endfor %}
              {% endfor %}

            "Score": "Score: {{ obj.Score }}"

          },
          "geometry": {
            "type": "Point",
            "coordinates": [{{ obj.latitude }}, {{ obj.longitude }}]
          }

        }{% if not loop.last %},{% endif %}
        {% endfor %}

How can I set and remove different coloured markers for certain scores this way?

I know about using jQuery to remove certain classes, but because I'm using jinja, I'm finding it difficult to understand how I would go about this.


texts: Extract lat long point values from EASE2 without reprojection

texts:
I have the following file:
EASE2
and want to extract the values for this geographic location (148.1, -35.6) #lon lat.
How can I do this without the need for reprojecting the ease grid file into geographic projection?

Here is an example code:

library(raster)
r <- raster (file)
e <- data.frame(148.1, -35.6)
ext <- extract(r, e)


texts: Raster mosaic not displaying raster values ArcGIS 10

texts:
I have eight DEM tiles that I downloaded from the USGS National Elevation Dataset.  The individual rasters are displaying correctly, but when I try to combine them in a Raster Mosaic I get arbitrary values (3.4+038 -> -3.4+038), which is considered to be NoData.

I've tried several times, half the time I get a warning message upon creation of the Mosaic (001003: Datum Conflict) even though I was sure to match coordinate systems.  I've even imported the coordinate system from one of the rasters going into the mosaic.

The mosaic shows the footprint of the DEM, but none of the data

I had a similar problem here, but I followed the process of converting to tif before I tried to mosaic the images.


texts: Iterating and extracting subdatasets from hdf files and exporting data into MS excel using ModelBuilder?

texts:
I have a large amount (10 years of data!) of MODIS LST .hdf files that I need to load into arc, extract subdatasets and import data into excel so it can be analysed. As I have such a large dataset I would like to automate this process as much as possibly.

I have limited experience of using ModelBuilder and python, but I'm willing to learn both. I have played around with ModelBuilder and from help from other Q&As (e.g. Using ModelBuilder to Iterate and extract subdatasets with in hdf file and project subdatasets?) and have built this:

Have been able to extract subdatasets and the pixels match up with the original image. I used this:

However, it seems that it runs the process for each .hdf file in the folder, but then it overwrites the output for each file so the final output is just one image (the last image processed).

Is there a way for the output to be a folder so that each image gets put into that folder?

I am now coming up with this error:

So it's saying the file extension does not support subdatasets, does anyone know how to solve this? I have tried all the file formats and none of them work!


texts: I'm doing quite heavy processing of Earth Observation data, and it appears that a lot of my workflow can be speeded up by using HDF5 for local storage. I mostly work with large stacks of raster files, and HDF5 as well as numpexpr allow me to very efficiently process large amounts of data. I would like to have as my external format HDF5, but my data can change geographic projection and geolocation. Ideally, I'd like GDAL to be able to use that information (not too worried about other tools), but I can't seem to find how to define that for these files.


texts: How to serve data from GeoServer on remote computer (server)?

texts:
I have developed an application using GeoExt + OpenLayers + GeoServer and it works fine. I have installed OpenGeo Suite in one of the computers at work (to serve as server) and I have uploaded my application online (cpanel).

The application runs well online but layers do not show because I am unable to connect to my server machine. I have tried to search for online resources but am rather confused on this issue.

I need help on this. Could someone guide me on how am supposed to set my server machine so that it can be accessible online so that layers and popups for my application can work?


texts: Rounding error in qgis point sampling tool

texts:
I'm not sure if this is a bug or of it is due to the implementation. I am using a raster that has values greater than 1,000,000 and Qgis displays these numbers in a scientific notation.  For instance, 1,101,702 gets displayed as 1.1017e+06 in the information dialogue.  I know that the raster stores this number correctly, because I polygonized it and the polygon got the correct identifier of 1101702.  But when I run the point sampling tool on the same area, the point gets marked as 1101700.  Does anyone know why this is?  And is there a simple way to fix this issue?

If not, can anyone recommend a way to script a simple point sampling tool that I can include in my data processing scripts?


texts: Calculating area of overlap between two polygons in separate layers using QGIS?

texts:
I have 2 layers with 1-1 polygon. The polygons are overlapping.

How do I calculate the area of the overlap using QGIS?


texts: How can i prevent a form from closing unless data is accepted by validation

texts:
I try to use a custom-form in QGIS 3 to validate input-data before saving the data in a database.
It would be nice to have the form open unless every input is correct.
But the dialog/form closes even if the code detects wrong data...

...
my_dialog = None
inputField = None

def form_open(dialog, layerid, featureid):
    global my_dialog
    my_dialog=dialog

    global inputField
    inputField = my_dialog.findChild(QLineEdit,"qle_input")

    buttonBox = my_dialog.findChild(QDialogButtonBox,"buttonBox")

    buttonBox.accepted.disconnect()

    buttonBox.accepted.connect(validate)

    buttonBox.rejected.connect(my_dialog.resetValues)

def validate():
    inputString=inputField.text() # get text from qleInput
    pattern = re.compile("[^A-Za-z0-9]"); #pattern of not allowed data
    if re.findall(pattern,inputString):
        #if match (wrong characters) - give a warning and let the dialog open to change text in qle_input
        textString='<b><center> Input contains wrong characters </center></b>'
        msgBox = QMessageBox()
        msgBox.setText(textString)
        msgBox.exec_()
    else: # Close dialog
        my_dialog.close


texts: Remove scene edge pixel error in annual composites in GEE

texts:
The median reducer over the cloud masked image collection is unable to remove the edge pixel error at a few scenes. Although the cloudmasking function has the code to remove such edge pixels not occurring in all scenes, but the median composites still has such errors at many location. The portion of composite with such error is attached along with the GEE code link for your kind reference and guidance.
I am getting such problems in a few more composites of another year. Please help me solve this issue.
https://code.earthengine.google.com/35fade3ad60b2634c1d9ec0d8fc748b5


texts:
I've been playing around with the GeoTools libraries, but I can only show shape files and no custom layers. Requirements are that it must use the SWTMapPane and not via the static SWTMapFrame.showFrame(). I'm looking for sample code that just creates features, say points, give it a style and show it on the pane.

Load shape file:

File openFile = new File(XXX);
if (openFile != null && openFile.exists()) {
    FileDataStore store = FileDataStoreFinder.getDataStore(openFile);
    SimpleFeatureSource s  =  store.getFeatureSource();
    Style style = Utils.createStyle(openFile, s);
    FeatureLayer featureLayer = new FeatureLayer(s, style);
    mapContent.addLayer(featureLayer);
    getMapPane().reset();
    getMapPane().redraw();
}

Now I want to add a custom layer containing features with a style, and then adding it to the mappane. Sample code I found always seem to contain deprecated code or JMapPane code.


texts: What is the Microstation keyin for Sychronize Google Earth View for a specific view?

texts:
I'm using TerraScan on top of Microstation V8i and I am setting up hotkeys for frequently used commands.  I can't find any documentation about the keyins for the Google Earth toolbar, specifically to link the Synchronize GE View button to a specific view (View 1).  I know the keyin for this tool is:

googleearth synch

I have followed some keyins with "selview1" or just the number 1 for View 1. Such as:

move up;selview 1
View Previous 1

Does anyone know the keyin to apply the synch google earth to a specific view?


texts: Do Personal Geodatabases show up when trying to Export Data using ArcGIS 10.1 SP 1 (on Windows 7)?

texts: I work with ArcGIS 10.1 SP 1 (on Windows 7), and I tried to export a feature layer to a personnal geodatabase using "right click in the table of content > data > export data". However, when I select "file and personnal geodatabase", only the file geodatabase are visible (though the mdb are visible in ArcCatalog). Of course, the obvious workaround is to use "feature class to feature class", but I would like to know why this "data > export data" is disabled. Is this because .mdb will disappear from ArcGIS in a near future ?


texts: I have a layer of line objects representing stream channels. These have been generated at 50K scale but do not have a Strahler order designation. The tutorials I've seen using the SAGA extension in QGIS all start at converting a DEM. I don't need this step. Is it possible to produce Strahler order on a vector layer of connected lines?


texts:
Using osm2po I found a problem. In the configuration file "osm2po.config" I can set the various values of speed limits for each road type but for example for this type of street (wtr.tag.highway.motorway = 1, 11, 120) in the output file the speed limit for that type of road is set to 110 instead of 120. Why is that?

I also wanted to know if openstreetmap provides data on speed limits for each road and if it is true, is there a way to insert them into the output file osm2po?


texts: How to handle NULL values in QgsExpression()?

texts:
When I try to add multiple conditions in QgsExpression()
eg:
QgsExpression(f"input='{abc}' OR value='{def}'") here abc = 1 and def = NULL. Even though abc is having value I am getting NULL in return. How to handle this?
Here abc and def are column names.
I need to get the row based on the value of abc even though def is NULL.


texts: qgis-web-client Not rendering label fonts

texts:
My install of qgis-web-client is not rendering label fonts, each character appears as a black box / square.

See screenshot:

Tested in FF & Chrome, have tried changing from the MS system font to Aerial with no change.

Thanks


texts: Exporting Feature Classes in SDE to shapefiles and tables

texts:
How can I export feature classes and tables in SDE to shapefiles and tables?

I am able to export the feature classes to shapefiles based on an extent. The problem I am running into is, I want to only export features from the tables that are related the one of the feature class. The tables and feature classes in the SDE are related based on an unique ID. If there are 30 features in the extent - I want to export only the related 30 records in the tables.

I am working in ModelBuilder.


texts: Plotting XY data from newly created shapefile in ArcMap?

texts:
I know how to plot XY data from a table and then convert to a shapefile, but here I'm trying to assign coordinates from a feature class's attribute table so that it would plot.

I created a point shapefile in ArcCatalog and used the Editor toolbar to enter some records manually. Two of these columns are for latitude and longitude.

How can I have ArcMap plot this layer using these coordinates?


texts: How to split name in label so only first two words appear?

texts:
The labels for my feature SSA_NAME are approximately 5 words long. I would like to shorten the labels to only 2 words. Every time I try a different code I receive either an error message or my layer is not labelled at all. Here is a code I tried based on another person's StackExchange question.

def convertLabel([SSA_NAME]):
         S= [SSA_NAME]
         S = S.split(" ")1
         return S


texts: Landsat8 surface reflectance rescaling pixel

texts:
Sorry for this basic question, but I am using Landsat8 images with "surface reflectance" treatment from Landsat8.

I am having some trouble with the data from some bands. Most of pixel are beyond 10k.

According Landsat8, we have a range between -2000 ~16000.
The valid range from Landsat8 is 0 ~ 10000, but I do not know what to do with the data outside this range

First option: -2000 will be 0
              16000 will be 10000
              And after this we multiply per 0.0001

Second option: lower than 0 will be -9999
               Higher than 10000 will be -9999
               And after this we multiply per 0.0001


texts: Merge and overwrite features between 2 shapefiles

texts:
I have 2 polygon shapefiles with similar features. One shapefile is "the master" (masterFile) and the other is a subset of features from the masterFile with some of those features slightly changed and/or a few new features added (subFile). I want to use the subsetFile to update the masterFile polygons.

What is the simplest way to do this?

I have already tried Merge and Spatial Join (perhaps I used them incorrectly), but the masterFile features that I wanted changed to the subsetFile features were not changed.


texts: Make a ground mask raster with the Greater Than tool

texts:
I'm trying to see if I can replicate the Esri solar suitability site example from this article.

I fail on the creation of the binary ground raster. I have made a Bare Earth raster from LiDAR ground points, and a surface raster from LiDAR first return points. The tool I use is Greater Than with first return as input raster 1, and ground raster as input raster 2.

I expect a raster with value 1 where the input1 raster is greater than the input2 raster, and value 0 where it is not greater. However, this is the result after the two input rasters:

The result should be a ground mask raster like this picture:

Okay, I have the raster I want after doing raster calc + 5 for the DTM and then using a CON expression concerning the DTM and DSM in raster calc to get the ground points raster.


texts: Speeding up Python Script tools get of raster properties using ArcPy?

texts:
I am using python script tools for ArcGIS 10.2, under Windows 8.1 64bit.

When I am using the arcpy.GetRasterProperties_management() command, it's very slow to process the data.

For example. I tried to get the maximum value for 14 rasters, if I use the batch mode of the system tool: Get Raster Properties, it uses only 0.19 seconds.

But when I am using python script tools to do this, it basically takes 1 second to print 1 value on the screen. So it takes 15 seconds to get the maximum value of the 14 rasters. So when I have hundreds of rasters to process, it takes too much time.

And this is when the CPU usage is 20%, which is the maximum usage when running ArcGIS on my laptop, sometimes the CPU couldn't even get there, when the usage is only 5 or 6%, it takes one minute to get all the values.

this is the code I use for this part:

for out_raster in out_rasters:
    try:
        raster_min = arcpy.GetRasterProperties_management(out_raster,"MINIMUM")
        raster_min = raster_min.getOutput(0)
        arcpy.AddMessage(raster_min)
    except:
        arcpy.AddMessage("No_data")

Is there any way I can increase the speed?