texts: NULL value not equal to empty value for attribute

texts:
I have a feature layer in a gdb. It has an attribute named uniqueid. Some of rows have a value for uniqueid and some do not. I am trying to select only the features that have a value (are not blank).

Using the Select by Attributes dialog box, I am trying to create an expression that will weed out the blank features. I tried using where: unique IS NOT NULL but that selected everything, both the features with a value and those that were blank. I then tried using Get Unique Values and chose '', which was one of the options. It through an error. What am I doing wrong here? How can I select just the features that have a value for the uniqueid attribute field?

Here's the actual text version of the error as requested:

  There was a problem selecting

  An invalid SQL statement was used. An invalid SQL statement was
  used.[fottprint_cslf_union_20181004] An invalid SQL statement was
  used.[SELECT OBJECTID FROM footprint_cslf_union_20181004 WHERE
  uniqueid IS NOT '']


texts: Algorithm to fill the gaps adjacent to parcels

texts:
I would like to find a way to fill the gaps between the road and its surrounding parcels, like described in picture below:

I have different parcels like A and B , i want to fill the gaps by adding the gap area to the adjacent parcel, the gap should be determined by the lines of the parcels which intersects the ROAD.

The final parcel B = (old) B + D , A = (old) A + C

I would like to find an algorithm to do this process.

I thought to make this algo :

Make a convex hull around ROADS and Parcels.
Erase Roads and parcels.
Get lines from the result layer and apply a tool like Extend line, but the constraint is that the lines would not be extended because they are already limited by parcel lines.

I'm not keen to follow this algorithm, i'm looking for any suggestions.

I'm using Arcpy in Arcgis 10.1, Arcobjects.


texts:
My QGIS 3.16 (and earlier 3.x versions) crashes when I open the "save as" dialogue box, even if the project is blank. It also seems really slow to run and load new data. QGIS 2.18.28 runs just fine on the same system: Dell Latitude laptop running Windows 10. I've uninstalled 3.x (but not 2.18) and reinstalled - it doesn't help. SAGA and GRASS bundled with 3.x seem to run ok, but haven't used them much. Any ideas?
Here's the Crash report:
Stack Trace

RtlpFreeHeapInternal :
RtlFreeHeap :
free :
DllUnregisterServer :
DllUnregisterServer :
CServerContextActivator::CreateInstance actvator.cxx:881
ActivationPropertiesIn::DelegateCreateInstance actprops.cxx:1983
CApartmentActivator::CreateInstance actvator.cxx:2189
CProcessActivator::CCICallback actvator.cxx:1640
CProcessActivator::AttemptActivation actvator.cxx:1527
CProcessActivator::ActivateByContext actvator.cxx:1371
CProcessActivator::CreateInstance actvator.cxx:1271
ActivationPropertiesIn::DelegateCreateInstance actprops.cxx:1983
CClientContextActivator::CreateInstance actvator.cxx:570
ActivationPropertiesIn::DelegateCreateInstance actprops.cxx:1983
ICoCreateInstanceEx objact.cxx:2045
CComActivator::DoCreateInstance immact.hxx:386
CoCreateInstance actapi.cxx:121
_SHCoCreateInstance :
CreateAndCacheIconOverlay :
CFSIconOverlayManager::GetOverlayIconInformation :
CFSIconOverlayManager::LoadIconOverlayIdentifiers :
CFSIconOverlayManager::_s_LoadIconOverlayIdentifiers :
CFSIconOverlayManager::LoadNonloadedOverlayIdentifiers :
EnableExternalOverlayIdentifiers :
CFSIconOverlayManager::GetFileOverlayInfo :
CFSFolder::_GetOverlayInfo :
CKnownFoldersFolder::_GetOverlayIndex :
CRegFolder::GetOverlayIndex :
CIconOverlayTask::DispatchTasks :
CIconOverlayTask::InternalResumeRT :
CRunnableTask::Run :
CShellTask::TT_Run :
CShellTaskThread::ThreadProc :
CShellTaskThread::s_ThreadProc :
ExecuteWorkItemThreadProc :
RtlpTpWorkCallback :
TppWorkerThread :
BaseThreadInitThunk :
RtlUserThreadStart :

QGIS Info
QGIS Version: 3.16.5-Hannover
QGIS code revision: 58ba7c1ed6
Compiled against Qt: 5.11.2
Running against Qt: 5.11.2
Compiled against GDAL: 3.1.4
Running against GDAL: 3.1.4

System Info
CPU Type: x86_64
Kernel Type: winnt
Kernel Version: 10.0.19042


texts: Failing to build PostGIS on OS X 10.7 Lion

texts:
I am trying to install postGIS from source to match the default install of postgreSQL on OS X 10.7. I installed GEOS and PORJ.4 from Kyngchaos and ran

$export MACOSX_DEPLOYMENT_TARGET=10.7
$export CFLAGS="-Os -arch i386 -arch x86_64"
$export PG_CPPFLAGS="-arch i386 -arch x86_64"
$export SHLIB_LINK="-arch i386 -arch x86_64"
$export LDFLAGS="-arch i386 -arch x86_64"
$cd ~/source/postgis
postgis $./configure --with-pgconfig=/usr/bin/pg_config --with-geosconfig=/Library/Frameworks/GEOS.framework/Versions/3/unix/bin/geos-config --with-projdir=/Library/Frameworks/PROJ.framework/Unix

I get this

PostGIS is now configured for x86_64-apple-darwin11.0.0

however after running make i get the following error

collect2: ld returned 1 exit status
lipo: can't open input file: /var/folders/dh/11gwx5cs5sl8v4r7qrzf3pd00000gn/T//ccsAdaB9.out (No such file or directory)
make[1]: *** [postgis-1.5.so] Error 1
make: *** [postgis] Error 2

Something i don't know about (architecture?) seems wrong here. Any clue?
Thanks in advance.


texts: Vertex editor in Qgis 3.0 : Snapping settings

texts:
Context :

I imported a grid point layer in Qgis 3.2 during my workflow, as seen in the picture below.

I now need to edit another point layer. As I don't need the first layer, I had it hidden, but I need to keep it in the project for possible future use.

Using the vertex tool editor, I want to move a point to specific coordinates. However, when I want to select the point in the second layer, all the points from the first layer are selectable, even if the layer is not editable. As the points are rather tight, it becomes hard to select the right point from the layer.

According to the research I made, this behaviour is not to be expected, as discussed in this link.

Question :

Is it a bug, or is there a setting that I am missing? I need to snap the vertex editor only to the active layer, any other behaviour feels really confusing.
Once the vertex is selected, is there another setting to enable or disable snapping to the hidden layers? In my mind, if a layer is hidden, it is meant to be forgotten until displayed again, so I would like to set Qgis according to that behaviour.


texts: How do you put GeoJSON attribute information in a popup with OpenLayers

texts:
I am having some trouble putting GeoJSON attribute information in a pop-up. The selection control is working, and I create an alert (for testing) and this is working. But I can't figure out how you can put GeoJSON attribute information in the pop-up.

The code I have so far is:

function init(){

        function onFeatureSelect(feature){
            alert("boom");
        }

        function onFeatureUnselect(feature){
            alert("unBoom");
        }

        var myStyles = new OpenLayers.StyleMap({
            "default": new OpenLayers.Style({
                fillColor: "${getColor}",
                strokeWidth: 1,
                strokeColor: "#000",
                fillOpacity: 0.8,
                graphicZIndex: 5,
                //label: "${getLabel}"
            },
            {
                context: {
                    getColor : function (feature) {
                        return feature.attributes.INW_T > 100000 ? '#ffc000' :
                               feature.attributes.INW_T > 10000 ? '#00317c' :
                                                                  '#FFEDA0' ;
                    }
                }
            })
        });

        var map = new OpenLayers.Map('map', {
            units: 'm',
            numZoomLevels: 19,
            controls: [
                new OpenLayers.Control.Navigation(),
                new OpenLayers.Control.ScaleLine(),
                new OpenLayers.Control.MousePosition(),
                new OpenLayers.Control.LayerSwitcher()
            ],
            projection: epsg900913,
            displayProjection: epsg4326 //Is used for displaying coordinates in appropriate CRS by MousePosition control
        });

        //LAYERS
        var lay_goo = new OpenLayers.Layer.Google('Google', {
            type: google.maps.MapTypeId.SATELLITE,
            sphericalMercator: true
        });

        var lay_osm = new OpenLayers.Layer.Google('OSM');

        var geojson_layer = new OpenLayers.Layer.Vector("GeoJSON", {
            styleMap: myStyles,
            projection: epsg4326,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "pc2.geojson",
                format: new OpenLayers.Format.GeoJSON()
            })
        });

        map.addLayers([lay_osm, lay_goo, geojson_layer]);

        selectControl = new OpenLayers.Control.SelectFeature(geojson_layer, {
            onSelect: onFeatureSelect,
            onUnselect: onFeatureUnselect
        });
        map.addControl(selectControl);
        selectControl.activate();

        map.setCenter(new OpenLayers.LonLat(lon, lat).transform(epsg4326, epsg900913), zoom);}

The GeoJSON has the following attributes:

{
"type": "FeatureCollection",

"features": [
{ "type": "Feature", "id": 0, "properties": { "PC2CODE": "10", "PC2NR": 10, "INW_T": 697830 }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 5.016855446988143, 52.349866914447176 ], [ 5.014431962156695, 52.347111770955848 ], [ 5.009350157428187, 52.344472149717674 ], [ 5.012133851358907, 52.344143065780877 ], [ 5.012314670330961, 52.343410242383534 ], more geometries

Hope someone can help!


texts:
I'm helping geocode addresses, and new at it. My office does not have Arc Server, is this required? I found this conversation Address or reverse geocode a point shapefile
that points me to address locators. Can I just load the street map and do this in ArcMap without ArcServer?

if not, I have been trying to do this through ArcGIS Explorer for Desktop, but it is suprisingly frustrating to migrate from this to ArcMap. It doesn't provide xy locations as a field, and only preserves one field from the original spreadsheet. In ArcMap I was able to do a join, and then 'add xy coordinates'. but there's got to be an easier way.


texts: Is there any way to know if a polygon is part of another bigger polygon(s) in SQL server?

texts:
Basicaly, what I am trying to do is we have a lot of geographical containers like territories, sub-territories etc. and the boundaries of each keep on changing all the times and I have to store the heirarchial data in to regular SQL using adjecent lists which says

  Territory A is parent of sub-territory B

Now, I was thinking that I can replace the entire thing if I can store each territory as a polygon in the map and make SQL responsible for finding which sub-territory is sub-geographical region of which Territory.

However, my knowledge of Geospatial data is next to basic and I wanted to know if SQL can do following for me:

Find out which territory(polygon) is a sub-territory(sub-polygon) of another one
If it can find which other territories have overlapping geographical definitions then that is fine too.


texts: Is there an advantage to installing an arcgis web adaptor in the DMZ

texts: I'm trying to architect a new ArcGIS Portal site. I see ESRI has instructions on configuring a reverse proxy to work with a web adapter installed on the network. However I have also seen sources talking about installing the web adapter directly on a DMZ web server and registering it with a portal on the network. Is there an advantage to doing it this way?


texts: Can I combine two models in ModelBuilder in ArcGIS?

texts:
I am using ModleBuilder within AcrGIS and am trying to create a modle that will:

1) Take my polygon layer and create an inverted buffer (using a negative value)
2) Add a field to the new polygon layer
3) Take that polygon layer and turn it into a raster
4) Take the raster and convert it to points
5) Add a field to the points layer
6) Use the spatial join tool to combine the attributes into one layer and table

I have made a model that encompasses all of the tools and analysis, but it has failed to run. It keeps reporting that it cannot complete the rater to point portion of the model. I have tried making two separate tools (one for the buffer, and another tool to create the points), and they ran without any problems. However that requires me to add the inverted polygon to the points model separately. It appears to only be an issue when the two are put together.

Am I doing something wrong, or is there a way for me to combine two separate models into one model?

The first image is the model with all of the components in it. The second image is of just the buffer model, and the third image is of creating the points.

Here is the code for the whole model I have built this far:

import arcpy

# Load required toolboxes
arcpy.ImportToolbox("Model Functions")

# Local variables:
B_mdb = "C:\\Users\\Desktop\\temp\\ModelBuilder\\B.mdb"
Name = "B"
B = "C:\\Users\\Desktop\\temp\\ModelBuilder\\B.mdb\\B"
inverted_buff_shp = "C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp"
inverted_buff_shp__2_ = inverted_buff_shp
bark_cells_tif = "C:\\Users\\Desktop\\temp\\ModelBuilder\\bark_cells.tif"
bark_points_shp = "C:\\Users\\Desktop\\temp\\ModelBuilder\\bark_points.shp"
bark_points_shp__2_ = bark_points_shp
inverted_buff_SpatialJoin = "C:\\Users\\Documents\\ArcGIS\\Default.gdb\\inverted_buff_SpatialJoin"

# Process: Iterate Feature Classes
arcpy.IterateFeatureClasses_mb(B_mdb, "", "", "NOT_RECURSIVE")

# Process: Buffer
arcpy.Buffer_analysis(B, inverted_buff_shp, "-0.35 Meters", "FULL", "ROUND", "NONE", "", "PLANAR")

# Process: Add Field (2)
arcpy.AddField_management(inverted_buff_shp, "model", "TEXT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Polygon to Raster
arcpy.PolygonToRaster_conversion(inverted_buff_shp, "OBJECTID", bark_cells_tif, "CELL_CENTER", "NONE", "2")

# Process: Raster to Point
arcpy.RasterToPoint_conversion(bark_cells_tif, bark_points_shp, "Value")

# Process: Add Field
arcpy.AddField_management(bark_points_shp, "orient", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Spatial Join
arcpy.SpatialJoin_analysis(inverted_buff_shp__2_, bark_points_shp__2_, inverted_buff_SpatialJoin, "JOIN_ONE_TO_ONE", "KEEP_ALL", "orient \"orient\" true true false 0 Double 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\bark_points.shp,orient,-1,-1;NO_NAME \"NO_NAME\" true true false 19 Double 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,NO_NAME,-1,-1;Shape_Leng \"Shape_Leng\" true true false 19 Double 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,Shape_Leng,-1,-1;Shape_Area \"Shape_Area\" true true false 19 Double 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,Shape_Area,-1,-1;BUFF_DIST \"BUFF_DIST\" true true false 19 Double 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,BUFF_DIST,-1,-1;ORIG_FID \"ORIG_FID\" true true false 10 Long 0 10 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,ORIG_FID,-1,-1;model \"model\" true true false 254 Text 0 0 ,First,#,C:\\Users\\Desktop\\temp\\ModelBuilder\\inverted_buff.shp,model,-1,-1", "INTERSECT", "", "")


texts: Show an esri leaflet feature's popup by default

texts:
I am trying to perform queries with esri leaflet and show the resulting polygon features with their appropriate popups by default (instead of having users click on the feature highlighted in order to open said popups). I have tried using the following lines, after the setting up the query parameters, with no luck:

map.addLayer(popup1); // where popup1 is a variable holding onEachFeature info

popup1.addTo(map);  // doesn't work either

map.openPopup(popup1);   // nothing

Note that the feature being queried (only one at at time) does show up fine, and clicking on it brings up the popups as expected. I just need these popups to be enabled by default, so users can close them if they want.


texts:
So I created my domain with coded values for my file geodatabase. I go to the feature class properties, fields tab, I create a new field, verify the data type of the field and domain matches and add that domain in the field properties below, click OK. Everything happens smoothly.

Now when I immediately reopen the same feature class properties I don't see the domain I set 1 minute ago. The domain field is empty and I can see in the ArcMap attribute table that the domain is not set to that field at all. And also now I can't edit or change any property of that field. It seems kind of locked.

This happens every time I set a domain to a field. Why is the domain disappearing? What am I doing wrong? The feature class is pre-created and I just want to add new field and add the domain to that field.


texts: leaflet / omnivore - create function to display omnivore CSV markers by property value

texts:
there seem to be a number of questions out there but none have quite answered this for me.

How would one go about taking a custom set of values from a property like latitude, and then assign them to change the marker color for omnivore.CSV data based on that value. Is there a way of calling the pointtoLayers functions within an onready style section and pointing it to fill colors based on these value tiers?

function getColor(d) {
    return d > 42.2 ? '#800026' :
           d > 39  ? '#BD0026' :
           d > 30  ? '#E31A1C' :
                      '#FFEDA0'; }

var MarkerOptions = {
    radius: 8,
    fillColor: getColor(feature.properties.latitude),
    color: "#000",
    weight: 1,
    opacity: 1,
    fillOpacity: 0.8 };

window.data3 = omnivore.csv('File.csv')
        .on('ready', function(layer) {
            this.eachLayer(function(markeroptions) {


texts: Openlayers with WMS service problem Resource interpreted as Image but transferred with MIME type application/vnd.ogc.wms_xml

texts:
I added a WMS layer

var wmslayer = new OpenLayers.Layer.WMS("Cities",{http://someip.com/ArcGIS/services/BurkeMcDowell2013_06_26FORGOPI/MapServer/WMSServer?request=GetCapabilities&service=WMS", { transparent: true },{isBaseLayer: false});}

And the URL which OpenLayers generated and call is:

{someip.com/ArcGIS/services/BurkeMcDowell2013_06_26FORGOPI/MapServer/WMSServer?request=GetCapabilities&service=WMS&SRS=EPSG%3A900913&TRANSPARENT=TRUE&VERSION=1.1.1&STYLES=&FORMAT=image%2Fpng&BBOX=-10018754.17,2504688.5425,-7514065.6275,5009377.085&WIDTH=256&HEIGHT=256}

But it's not returning images but a message in console:

Resource interpreted as Image but transferred with MIME type application/vnd.ogc.wms_xml: http://someip.com/ArcGIS/services/BurkeMcDowell2013_06_26FORGOPI/MapSe
X=-10018754.17,009377.085,-7514065.6275,7514065.6275&WIDTH=256&HEIGHT=256". OpenLayers.js:460)
OpenLayers.Tile.Image.OpenLayers.Class.setImgSrc OpenLayers.js:460
OpenLayers.Tile.Image.OpenLayers.Class.initImage OpenLayers.js:460
OpenLayers.Tile.Image.OpenLayers.Class.renderTile OpenLayers.js:456
OpenLayers.Tile.Image.OpenLayers.Class.draw OpenLayers.js:456
OpenLayers.Layer.Grid.OpenLayers.Class.initGriddedTiles
OpenLayers.js:480
OpenLayers.Layer.Grid.OpenLayers.Class.moveTo OpenLayers.js:469
OpenLayers.Layer.OpenLayers.Class.redraw OpenLayers.js:265
OpenLayers.Layer.HTTPRequest.OpenLayers.Class.redraw
OpenLayers.js:451
OpenLayers.Map.OpenLayers.Class.addLayer OpenLayers.js:230
Objs.adddlayer CreateMapMediator.js:1531
f.event.dispatch jquery-1.7.2.min.js:3
h.handle.i

Can anyone help me solve it? (I added someip.com coz I can't add a numeric IP there)


texts: I'm curious to know what is the counterpart for ArcGIS focal statistics tool in QGIS?


texts: Raster appears smooth when created but pixelated when opened again

texts:
I have a set of xyz points that I want to turn into a raster format for bathymetry.  I am using the Topo to Raster (3D analyst) tool and using the depth as the Field and PointElevation as they Type.  When the tool finishes it outputs a very smooth looking raster.  However when I open a new ArcMap document and pull in the created raster it looks very pixelated.  I can't figure out what is causing the change.  One thing I did notice is that at the bottom of the original Arcmap window the units are listed in Decimal Degrees while in the new Arcmap the units are listed as meters.

Can anyone explain to me what is happening when I am opening the raster in a new ArcMap window?  How can a switch the dataframe between meters and decimal degrees and will this change fix the pixelation?


texts: Trimble GeoXH, XM, XT - Data logging without Terrasync or GPSCorrect?

texts: Is it possible to capture data for post processing on a Trimble GeoXH, XT, or XM without using Terrasync or GPScorrect? Specifically I would like to create a RINEX file. A Trimble 4000 series receiver can log data internally, that can then be exported to a PC and converted into a RINEX file. Is this possible to do with a stock Geo XT, XH, XM?


texts: Is there a shortcut for switching between drawing Lines and drawing Arcs in ArcMap?

texts:
I normally use AutoCAD for something like this.

I am digitizing parking lots and have to press the buttons on the toolbar to switch between lines and arcs while drawing.

Is there a shortcut that I can type rather than clicking the buttons?


texts: Where to find and share Python geoprocessing scripts?

texts: Where are folks storing and sharing their Python geoprocessing scripts?  For example, I know that GitHub has some shared scripts.  ESRI is phasing-out ArcScripts, which is still really handy, and directing people to the Geoprocessing Model and Script Tool Gallery.  What other repositories are available for sharing geoprocessing scripts?


texts: Values within an attribute table are displayed in a legend in a layout using PyQGIS

texts:
I currently have a script that runs data within my map and then exports it as a layout in print composer. On the map, I have a legend key that is reference to the data on the map. Within my legend, I have a yellow buffer that is titled 'MM1 1.0 Buffer'. Within the attribute table of this vector file, I have a column called 'LAC_Population_sum', with a value of affected population from an earthquake, within that buffer area. However, I can not seem to find a way to get this value, so it is displayed within my legend on the layout.
The script I use to create a memory buffer is:
# Creating the distance buffers
layer = iface.activeLayer()
feats = [ feat for feat in layer.getFeatures() ]
epsg = layer.crs().postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer&field=x:real&    field=y:real&field=point_id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                       'MMI 1.0 Buffer',
                       'memory')
prov = mem_layer.dataProvider()

#Adding distance to the buffer
for i, feat in enumerate(feats):
    point = feat.geometry().asPoint()
    new_feat = QgsFeature()
    new_feat.setAttributes([i, point[0], point[1], feat.id()])
    tmp_feat = feat.geometry().buffer(DistInt[1.0], -1) #printing a MMI 1.O Intensity Buffer for a distancethat I have calculated previously and is stored within a dictionary (DistInt)
    new_feat.setGeometry(tmp_feat)
    prov.addFeatures([new_feat])

QgsProject.instance().addMapLayer(mem_layer)

My Attribute table for the MM1 1.0 Buffer, with the value of LAC_Population_sum:

What my legend currently shows:

My code that creates the LAC_Population_Sum field:
#Using zonal statistics to work out affected population sum column
root = QgsProject.instance().layerTreeRoot()
for layer in root.children():
  if layer.name().startswith('LAC'):
    prefix = layer.name()
    params = {'INPUT_RASTER': layer.name(), 'RASTER_BAND': 1, 'INPUT_VECTOR': 'MMI 1.0 Buffer', 'COLUMN_PREFIX': prefix+'_', 'STATS': 1}
    processing.run("qgis:zonalstatistics", params)


texts: Feature classes not visible to some ArcSDE users?

texts:
Some of my organization's data (several SDE databases) has been migrated to a new server and we have recently upgraded from ArcGIS 10.2 for Server to ArcGIS 10.4 for Server and SQL Server 2005 to SQL Server 2014. Users are using ArcGIS for Desktop 10.1 or 10.2.

Now when some users attempt to view the migrated data, no feature classes are shown. Feature datasets are shown however. For others, the feature classes are present and may be worked with as usual. Any ideas of why this is happening and what may be done to allow all users to see these feature classes?

Feature classes have been absent for some users using Desktop versions 10.1 and 10.2. So I do not think it is a Desktop version issue. I do not think it is a Native Client version issue because I believe users would be prevented from viewing the database entirely.


texts:
I have a problem with hillshade function: I've set 315 azimuth, 45 vertical angle and 1 Z factor, but the result is very bad...
The original dem file is a .geotiff from Global Data Explorer (Aster Global DEM V2). I would like to obtain a pretty hillshade, but I don't know what is wrong.

Here the original and the screen of the hillshade result. It's 1:660.000 scale


texts: Using Update cursor from data access module?

texts:
I am currently attempting to use an update cursor from the data access module to update the PPCVal field based on a certain a set of conditional logic rules that I've laid out below.  I've included the structure of my table  so you can see the field names:

[![enter image description here][1]][1]

I've pasted a screenshot of my code below:

[![enter image description here][2]][2]

The code executes without any exceptions or errors but it appears it is not executing correctly on my second "elif" conditional where features that have WATER_FLAG = Z where index 3 (PPCVal) should be equal to index 1 (PPC_LOW).

You can see in my first screenshot of the table that the PPCVal field (index 3) is not updating correctly with the PPC_LOW value when WATER_FLAG = Z.

import arcpy, os

arcpy.env.workspace = r'C:/Capstone/Vermont1.gdb'
arcpy.env.overwriteOutput = True

fc = 'C:/Capstone/Vermont1.gdb/VT_PTS2_Intersect3'
fields = ['WATER', 'PPC_LOW', 'PPC_HIGH', 'PPCVal', 'WATER_FLAG']

with arcpy.da.UpdateCursor(fc, fields) as cursor:
    for row in cursor:
    #Evaluate WATER field for 'Y' and apply value from PPC_LOW (index 1)
    #to PPCVal
    if (row[0] == 'Y'):
        row[3] = row[1]
    #Evaluate WATER field for blank ' '  and apply value from PPC_HIGH
    # (index 2) to PPCVal
    elif (row[0] == ' '):
        row[3] == row[2]
    #Evaluate WATER_FLAG field for 'Z' (hauled water) and apply value from
    # PPC_LOW (index 1) to PPCVal
    elif (row[4] == 'Z'):
        row[3] = row[1]
    #Otherwise PPCVal = PPC_High (index 2)
    else:
        row[3] = row[2]

    cursor.updateRow(row)

print "PPCVal update finished."


texts: Stitching thousand of 32 bit Tiff without geographic information

texts:
I'm working on a project which goes beyond my knowledges, I need your help.

I'm working with FLIR movies (thermal imagery recorded from plane), and I have to find a way to make a Tiff containing the temperature information and geographic information.
I can extract images from the movies in different ways (png, jpeg, csv, Tiff, 32 bit tiff, FITS).

My plan was to extract one jpeg/second and then stitch those thousand of pictures with hugin, then georeference the stitched result on ArcGis.

My problem is that I suppose the JPEG will be in RGB format so the temperature won't be displayed...

So the only way might be to export 1 32 bit tiff/second (which is the only format with the real temperature information) and then try to stitch them together... But I don't see how.
The study area is really huge so I have 5000 thousand pictures each day, so I can't afford any manual method.

Does someone has a genuis idea?


texts: Retrieve all files created by a task to my local folder

texts:
I'm exporting a map to my drive using the GEE python API, I'm looking for a way to retreive all the created files automatically to a local folder.
my current code create a image on my drive :
def download_to_disk(country_code, year):
    """download the tile to the GEE disk"""
    file_name = 'alerts_' + country_code + '_' + year
    task_config = {
        'image':alerts,
        'description':file_name,
        'scale': 30,
        'region':country,
        'maxPixels': 1e10
    }

    task = ee.batch.Export.image.toDrive(**task_config)
    task.start()

Is there a way to retrieve the files created by this very task automatically using only the Python API ?


texts: Can't reseed tiles on mapproxy from geoserver wms

texts:
I am confused on how to initiate reseeding in mapproxy. I used the file mapproxy.yaml and seed.yaml for the first seeding. Its source has been significantly updated. So I want to cache everything anew. This mapproxy-seed -f mapproxy.yaml -s seed.yaml doesn't work (tiles weren't created like  with the dry run setting). I failed with this too:

--reseed-interval 90d --reseed-file reseed.time \
--continue --progress-file .mapproxy_seed_progress

How can I reseed tiles?
This my yaml file:

# ---------------------------------------
# MapProxy example seeding configuration.
# ---------------------------------------
#
# This is a minimal MapProxy seeding configuration.
# See full_seed_example.yaml and the documentation for more options.
#

seeds:
  myseed1:
    caches: [osm_cache]
    # grids: []
    # coverages: []
    levels:
      to: 12
   # refresh_before:
    #  time: 2013-10-10T12:35:00

#cleanups:
 # myclean1:
  #  caches: [osm_cache]
   # remove_before:
    #  days: 14
    #levels:
     #   from: 11

#coverages:

# -------------------------------
# MapProxy example configuration.
# -------------------------------
#
# This is a minimal MapProxy configuration.
# See full_example.yaml and the documentation for more options.
#

# Starts the following services:
# Demo:
#     http://localhost:8080/demo
# WMS:
#     capabilities: http://localhost:8080/service?REQUEST=GetCapabilities
# WMTS:
#     capabilities: http://localhost:8080/wmts/1.0.0/WMTSCapabilities.xml
#     first tile: http://localhost:8080/wmts/osm/webmercator/0/0/0.png
# Tile service (compatible with OSM/etc.)
#     first tile: http://localhost:8080/tiles/osm/webmercator/0/0/0.png
# TMS:
#     note: TMS is not compatible with OSM/Google Maps/etc.
#     fist tile: http://localhost:8080/tms/1.0.0/osm/webmercator/0/0/0.png
# KML:
#     initial doc: http://localhost:8080/kml/osm/webmercator
services:
  demo:
  tms:
    use_grid_names: true
    # origin for /tiles service
    origin: 'nw'
  kml:
      use_grid_names: true
  wmts:
    # use restful access to WMTS
    restful: true
    # this is the default template for MapProxy
    restful_template: '/{Layer}/{TileMatrixSet}/{TileMatrix}/{TileCol}/{TileRow}.{Format}'
    # and also allow KVP requests
    kvp: true
  wms:
    md:
      title: MapProxy WMS Proxy
      abstract: This is a minimal MapProxy example.

layers:
  - name: osm_background
    title: osm_background
    sources: [osm_cache]

caches:
  #geoserver cache:
  gasps_osm_cache:
    # cache the results in two grids/projections
    grids: [osm_grid]
    sources: [osm_source]

sources:
  gasps_osm_source:
    type: wms
    req:
      url: 'http:/xxx.local:xxxx/geoserver/wspace/wms'
      layers: wspace:layer
grids:
    osm_grid:
        base: GLOBAL_MERCATOR
        srs: 'EPSG:3857'
        origin: nw

globals:


texts: I can't seem to find an answer with regards to obtaining COGO for ArcGIS Basic 10.2. Does COGO fall under any of the extensions currently offered by ESRI or does it only fall under the higher level licenses.


texts:
I have this list:

shapefiles that has these directories for the files: (pairs of pst and dbound) in each folder.

[['C:\\Users\\user\\Desktop\\eff\\20194\\DBOUND\\DBOUND.shp',
  'C:\\Users\\user\\Desktop\\eff\\20194\\PST\\PST.shp'],
 ['C:\\Users\\user\\Desktop\\eff\\20042\\DBOUND\\DBOUND.shp',
  'C:\\Users\\user\\Desktop\\eff\\20042\\PST\\PST.shp'],
 ['C:\\Users\\user\\Desktop\\eff\\20161\\DBOUND\\DBOUND.shp',
  'C:\\Users\\user\\Desktop\\eff\\20161\\PST\\PST.shp'],
 ['C:\\Users\\user\\Desktop\\eff\\20029\\DBOUND\\DBOUND.shp',
  'C:\\Users\\user\\Desktop\\eff\\20029\\PST\\PST.shp'],
 ['C:\\Users\\user\\Desktop\\eff\\20008\\DBOUND\\DBOUND.shp',
  'C:\\Users\\user\\Desktop\\eff\\20008\\PST\\PST.shp']]

When I try this to export the shp returns:

  ValueError: Geometry column cannot contain mutiple geometry types
  when writing to file.

import geopandas
for i in shapefiles:
    pst = gpd.read_file(i[0])
    dbound = gpd.read_file(i[1])
    dbound.reset_index(inplace=True)
    wdp = gpd.sjoin(pst, dbound, how="inner", op='within')
    wdp['DEC_ID']=wdp['index']

    folder=i[0].split("\\")[-3]
    o = r'C:\Users\user\Desktop\yyyyyh\wdp_'+folder
    wdp.to_file(o)

but when I try it manually (without the loop) it works perfect. So how these geometry types are a problem only with the loop. I literally tried the same shapefiles and all fine.


texts: Repeating labels even if I set to "No repeat" in Placement

texts:
In QGIS, I'm working on a map with the contours of a certain location. For the contour labels I used a fixed position with the "Label Toolbar". In the settings of the labels, in "Placement" I set it to "No repeat". Even so some labels are still repeating, in Print Layout as well. Now how can I place only one label for each feature in my drawing?
Repeating labels in the project:

Repeating labels in Print Layout:


texts: In ArcMap, how to create random points of a specific amount in each polygon?

texts:
I am doing an accuracy assessment of land-use product and need random points as test points. What I'm trying to do is to generate different numbers of random points in different polygons (like 5 for smaller ones and 10 for bigger ones). But in ArcMap, the tool to do this can only generate a specific total amount of points, or to create points according to a specific field for the constrained feature class.

I tried to create a field called 'RP amount' and use field calculator to do some math about the amount of random point using area data, but there are many polygons belong to one parcel, and it is difficult for those polygons to have an evenly distributed random points.


texts:
Can I use a greenest pixel or a no clouds Sentinel 2 for NDVI analysis? I want to use it for extract spectral signature.
And now that I´m asking about it, I want to know if I can work classification based on spectral signature in Google Earth Engine.

Note: Sorry if I have a bad redaction, my native language is note english.


texts: Is it possible to return coordinates as arrays rather than strings using ogr?

texts:
I am trying to convert a list of coordinates between CRS's using gdal's Python bindings. I would like to have the converted coordinates output as an array, but I can only seem to export as various types of strings.

import ogr
import osr
import numpy as np

xs = np.array([11705274.6374, 10705274.6374])
ys = np.array([4826473.6922, 4856473.6922])

# Spatial Reference System
inputEPSG = 3857
outputEPSG = 4326

# Create an ogr object of multipoints
points = ogr.Geometry(ogr.wkbMultiPoint)
for i in range(len(xs)):
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(float(xs[i]), float(ys[i]))
    points.AddGeometry(point)

# Create coordinate transformation
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(inputEPSG)

outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(outputEPSG)

coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# transform point
points.Transform(coordTransform)

points.FlattenTo2D()
print(points.ExportToWkt())

MULTIPOINT (105.15027111593 39.7278572772792,96.1671182747348 39.9348108861792)

The xs and ys provided are just dummy points I randomly typed. I am casting them as numpy arrays as my program will be providing them as such.

I can break the MULTIPOINT string up and pull the coordinates, but this is undesirable as it may cost precision and is just clunky. Any ideas?


texts: Processing many NetCDF files using ArcPy gives Serious Application Error

texts:
I am processing more than 5k NetCDF files in ArcPy, using ArcGIS Desktop 10.5. I am doing it using the Python console of ArcMap.
The process runs well, but when around 113 files are done, the process breaks down. The ArcMap gives an ERROR(attached at the end).
Does anyone know if there is a limitation with the number of files to process?
The process is:

Read the netCDF files --> arcpy.MakeNetCDFRasterLayer_md

Apply Zonal Statistics as table -->arcpy.gp.ZonalStatisticsAsTable_sa

Export the results to a CSV

The code is:
in_fd="is the netCDF files directory. More than 5k files."
shape= "is the shapefile where the zonal statistics will be applied."

names=[]
path=[]
dic={}

for root,folder,files in os.walk(in_fd):
for file in files:
    if os.path.isfile(os.path.join(root,file)) and file.endswith('.nc'):
        path.append(os.path.join(root,file))
        names.append(os.path.basename(file))

arcpy.env.workspace = r"Z:\EMSV-064_IM_for_EMSN071_ForestFireHazard_Germany\5_Process\3_Task2\1_Data\4_FWI\MIKEL\FWI_mv.gdb"
arcpy.env.overwriteOutput = True
for i in range (0,len(path)):
    print("#############")
    print("Empieza el:   " + names[i])
    ras=arcpy.MakeNetCDFRasterLayer_md(in_netCDF_file=path[i], variable="dc", x_dimension="longitude", y_dimension="latitude", out_raster_layer=names[i][:-3], band_dimension="", dimension_values="", value_selection_method="BY_VALUE")

    table=arcpy.gp.ZonalStatisticsAsTable_sa(shape, "Id", ras, os.path.join(r"Z:\EMSV-064_IM_for_EMSN071_ForestFireHazard_Germany\5_Process\3_Task2\1_Data\4_FWI\MIKEL\Statistics_DC","id_"+str(names[i][13:21])), "DATA", "ALL")

    with arcpy.da.SearchCursor(table,"*") as cursor:
        for row in cursor:
            dic[names[i][13:21]]=row[4],row[5],row[7],row[8]

    ls=[]
    vls=[(k,)+v for k,v in dic.items()]
    for row in vls:
        ls.append(row)

    if os.path.exists(r"Z:\EMSV-064_IM_for_EMSN071_ForestFireHazard_Germany\5_Process\3_Task2\1_Data\4_FWI\MIKEL\CSV_mv\FWI_DC.csv"):
        os.remove(r"Z:\EMSV-064_IM_for_EMSN071_ForestFireHazard_Germany\5_Process\3_Task2\1_Data\4_FWI\MIKEL\CSV_mv\FWI_DC.csv")
    with open(r"Z:\EMSV-064_IM_for_EMSN071_ForestFireHazard_Germany\5_Process\3_Task2\1_Data\4_FWI\MIKEL\CSV_mv\FWI_DC.csv","wb")as out:
        csv_out=csv.writer(out)
        csv_out.writerow(['Fecha','MIN','MAX','MEAN','STD'])
        for row in ls:
            csv_out.writerow(row)

    arcpy.Delete_management(ras)
    arcpy.Delete_management(table)

The error is:


texts: Client side caching in ArcGIS JSAPI maps

texts:
Is it possible to cache a map layer or feature layer on client side (the browser) using ArcGIS JSAPI, if so, can someone provide a link to documentation and/or samples?

If not, what other client technologies (e.g. Flex, .Net, SilverLight) support caching in ArcGIS?


texts: How to not overwrite with NoData values when updating a raster file using a numpy array?

texts:
I have a numpy array (called out_window) with the following values:

out_window:

[[[499 498 481]
  [494 490 464]
  [483 480 -1]
  [482 -1 -1]]]

This numpy array was obtained by clipping a raster file using a polygon feature and altering only those elements with value different to -1 (-1 is NoData)

I'm trying to update that same raster file with those new values.

import rasterio
import geopandas

raster_in = r'C:\.......\raster_in_file.tif'
raster_out = r'C:\.......\raster_out_file.tif'
copy(raster_in, raster_out)
polyg = gpd.read_file(gdb_vml_polygs, driver='FileGDB', layer=fc_vml_polygs)

out_window, out_transform = mask(raster_in, polyg['geometry'], all_touched=True, crop=True)

out_meta = raster_in.meta.copy()

#this is not relevant
width = ...
height = ...
x = ...
y = ...

with rasterio.open(raster_out, "r+", **out_meta) as dest:
    dest.write(out_window, window=Window(x, y, width, height))
dest.close()

The problem is that the instruction dest.write(.....) is overwriting the raster with all 12 values of the out_window array. I do not want to overwrite the raster with the -1 values. Those 3 elements of the raster should be left unaltered.

Is there any way to do this?


texts:
I have a FeatureLayer that uses a ClassBreaksRenderer. The renderer loads data from a collection like so:

var renderer = new esri.ClassBreaksRenderer(
            null,
            function (graphic) {
                return data[graphic.attributes.ID];
            });

Using a LabelLayer, is there a way to only render labels when the feature has a value (as determined by the ClassBreaksRenderer)?

I am able to "hide" the labels in this instance, by setting the font size to 0, but the DOM elements are still created. On a layer with thousands of features where only a sub-set are displayed, that kills performance.


texts: Restrict geoprocessing from ArcGIS Server service

texts:
I have an Image Service up on our ArcGIS server. Its a DEM we are providing it to some of our partners. Works great and fast, but I want to turn off the ability to run tools from the service within ArcGIS programs.

I never knew that you could, but tried creating contours, exporting, it all works and I need to restrict that. The reasoning is because anyone can create copies and derivatives of this information and we don't want that. I've tried removing mensuration and metadata in the capabilities, but it doesn't appear to do anything.

Any ideas? Google has not helped me out yet.


texts: Measure each angle of a selected polygon

texts:
This post is also a continuation of a previous one. I have the code (in VBA and ARcObjects) shown below to measure each angle of a selected polygon in degrees.
I have the following problems:
1. Some angles are negative and I do not wish this since I measure the interior angles of a polygon.
2. Some angles are measured wrong i.e. I get the reflex angle instead of the normal angle. e.g. 270 instead of 90 degrees.

Could somebody see the code and help me for the neccessary modifications please?
There is a procedure that calls a function shown also below. Thanks.

    'This procedure gives the angles between the vertices of a polygon

' ----------------------------------------

Dim pMxDoc As IMxDocument
Dim pGeoFtrLyr As IFeatureLayer
Dim pFtrCls As IFeatureClass
Dim pFtrCsr As IFeatureCursor
Dim pFtr As IFeature
Dim pGeomColl As IGeometryCollection
Dim pSegColl As ISegmentCollection
Dim pCurve As ICurve
Dim pLine1 As ILine
Dim pLine2 As ILine
Dim dVtxAng As Double
Dim i As Long, j As Long
Dim PI As Double

'    ' Calc PI
PI = Atn(1) * 4

' Get a cursor on the features in the first layer
Set pMxDoc = ThisDocument
Set pGeoFtrLyr = pMxDoc.FocusMap.Layer(12)
Set pFtrCsr = pGeoFtrLyr.Search(Nothing, False)

'Select a polygon

Dim pEnumFeature As IEnumFeature
Set pEnumFeature = pMxDoc.FocusMap.FeatureSelection
pEnumFeature.Reset
Set pFtr = pEnumFeature.Next
If pFtr Is Nothing Then
   MsgBox "Please select a valid feature", vbInformation
 Exit Sub
End If
If Not TypeOf pFtr.ShapeCopy Is IPointCollection Then
MsgBox "Not a valid feature", vbInformation
Exit Sub
End If

    ' Get all disjoint paths
    Set pGeomColl = pFtr.ShapeCopy

    ' Loop thru the paths
    For i = 0 To pGeomColl.GeometryCount - 1
        ' Get all segments that make up this path
        Set pSegColl = pGeomColl.Geometry(i)
        Set pCurve = pSegColl

        ' Loop thru the segments
        For j = 0 To pSegColl.SegmentCount - 2

            ' Get the next two lines
            Set pLine1 = pSegColl.Segment(j)
            Set pLine2 = pSegColl.Segment(j + 1)

            ' Calculate the left side angle (in degrees)
            dVtxAng = CalcAngleBetweenLines(pLine1, pLine2, True)
             MsgBox Abs(dVtxAng * 57.2957795)

        Next j

        ' Check for a closed polylin

        If pCurve.IsClosed Then
            Debug.Print pFtr.OID; " is closed!"
            Set pLine1 = pSegColl.Segment(j)
            Set pLine2 = pSegColl.Segment(0)
            dVtxAng = CalcAngleBetweenLines(pLine1, pLine2, True)

            'get each angle of the polygon in degrees
            MsgBox dVtxAng * 57.2957795

        End If

    Next i

   End Sub

Private Function CalcAngleBetweenLines(pLine1 As ILine, pLine2 As ILine, AsDegrees as Boolean) As Double

Dim dAng1 As Double, dAng2 As Double, dVtxAng As Double, PI As Double

dVtxAng = pLine1.Angle - pLine2.Angle

' Get the reversed angle of the first line
pLine1.ReverseOrientation
dAng1 = pLine1.Angle
If dAng1 < 0 Then dAng1 = 2 * PI + dAng1
pLine1.ReverseOrientation

' Get the angle fo the second line
dAng2 = pLine2.Angle
If dAng2 < 0 Then dAng2 = 2 * PI + dAng2

' subtract 2 from 1 - that's the smallest angle between the two lines
dVtxAng = dAng1 - dAng2

' Ensure we have the angle to the left of the lines

If dVtxAng >= 2 * PI Then
dVtxAng = dVtxAng - 2 * PI
ElseIf dVtxAng < 0 Then
dVtxAng = dVtxAng + 2 * PI
End If

CalcAngleBetweenLines = dVtxAng
End  Function


texts:
I am using the clipper tool to cut a raster layer to a shapefile. I have done this successfully to 2 raster layers. Now I am working on a third raster layer from a different source and when I hit okay to clip the layer, it just acts like it is trying to load forever. It never completes the process, and it never throws an error.

Why would this be?

For the layer that clipped correctly, the CRS is w6884, and for the one that didn't work, it is nad83

It seems that any operation I try to perform on this raster layer never completes, yet doesn't throw an error. Is there something wrong about the layer I am using?

The file that works was .bil and the one that didn't work is a .adf. I'm not sure what the difference is.


texts:
I have a set of polygons in ArcGIS 10, which I want to label. Despite setting the Unplaced label icon to "on" and using an expression for the label (so that all objects should be labelled), not all objects are being labelled, whether using the standard label engine or Maplex.

I have tried reloading Arc, creating annotations, creating a new feature class, but nothing works. Has anyone else seen this behaviour?

There is a similar question How to "Never Remove" labels in Maplex? but this seems to be a slightly different issue - I can't even get the labels to show up as unplaced.

The screenshot below gives some more context - I would expect all the yellow regions to have a black or red (denoting unplaced) A, but they don't.


texts: ArcGIS Server 10.1 - Inconsistent Querying Errors

texts:
I have recently upgraded to ArcGIS Server 10.1 SP1.  Since making the switch, I am getting a lot of inconsistent errors, specifically performing queries against a map service's REST API with the Javascript API.

I say inconsistent because an identical query will sometimes work, and other times I get an 'unable to complete operation' error.  The queries are extremely simple, for example the where clause is "COUNTRY_ID = 2", and I am positive its not the where clause because it always worked coming from ArcGIS Server 10.0.

Testing the query in the REST page directly from the services directory yields the same result.  My result alternates between success and error - query 1: success, query 2: error, query 3: success, query 4: error.

Has anyone experienced something similar?

UPDATE:  From info I found in this post, I found that making POST requests resolves my issues.


texts: I would like my map to show more color differences between sites that have 1-200 turbines. I have 1-293 turbines per location. However I would like to show the areas that have higher numbers as having more of an impact. For example the area that is in the upper right hand side is a very clustered area and a bird is more likely to hit a turbine in an area of high density. I hope that makes sense. Sorry if I need to explain more let me know. I am using ArcGIS 10.1 and my data came from USGS wind energy website. The attributes within that I am interested in are location, and total number of turbines. However my problem is that with such a large area I am having a large amount of over lap and unneeded noise. I am trying to figure out the best way to show large farms or clusters of farms without 20 points all within the same location.


texts: PostGIS - Interpolate point between contours

texts:
Using PostGIS, what is the best method to interpolate the z value for points (the orange dot in the image below) located between each contour layer (blue lines and red lines in the image below)? So for the blue contours it's going to be ~ 917 and for the red contours it'll be ~ 819.

The blue contours are a 2D Multipolygon and the red contours are a 2D Multilinestring (so for both, the contour value is an attribute field. If necessary I could update the geometries to have a z value).

Thanks!


texts: Selecting all features except one using QGIS?

texts:
Using QGIS 2.18. I have 346 polygons on a layer and I want to display all but one.

In this instance I wish to run a query using field ONSCode but show all except 17UB.

I have been trying to use NOT but in the wrong context and get errors.


texts: Map of Spratley and Scarborough Shoal EEZs

texts:
I have scoured the internet looking for a map of the South China Sea that shows the disputed Spratly Islands and Scarborough Shoal along with the EEZ (200 nm) associated with them.
In other words, one of the big disputes of course is whether these features are "islands" which give the 200 nm exclusive economic zones (EEZ) all around. Assuming that the features are islands (and not low-tide elevations or rocks), I'm looking for a map that shows the islands plus their putative EEZ.
You can find lots of maps of the South China Sea, and they will often include markings for the EEZs of the countries. But the hot dispute is over a) who owns the features in the ocean, and b) whether you can get an EEZ go to with it. I'm interested in the latter question and I'm looking for a map to show what they'd look like if indeed recognized.


texts:
I'm looking for a universal world grid of 5X5 km or 10X10 kms.
Is there any known reference in place? And is it possible to get the data ?


texts: Select points on line by coordinates ModelBuilder?

texts: I would like to use ModelBuilder to select points along a polyline between two sets of coordinates then export it as a shapefile (the exporting I can do). Each points represents 1 ft and the line represents a stream. This would be used to calculate sinuous length of stream no flows or areas of interest. So pretty much the number of points between X1,Y1 and X2, Y2 would be the length of the no flow or AOI.


texts: Smart way of finding the polygon intersection of a huge GeoJSON file with a small GeoJSON file without looping through all polygons

texts:
I have a huge GeoJSON file (over 300MB) which has multiple different polygons. Also, I have a very small GeoJSON file with only one polygon that doesn't have more than ten points. Using the intersection function from turf.js, I can create a simple logic where I loop through all polygons inside my huge GeoJSON file and I try to find which of its polygons intersect with my small GeoJSON file. This is a functional logic, but I'm not sure if that's how I'm supposed to solve this kind of problem.
What would be a smarter way of solving this problem?
How can I find which of my polygons (inside a huge file) intersect with a small polygon (inside my small file) without searching if there's an intersection in all polygons?
Even though I'm using turf.js with node.js I'm not stuck with this library and language. So, if there's a better way of solving it with python or command line I can change my approach for solving this problem. I think I'm missing a concept for solving problems with this kind of scale.


texts: ArcPy script not finding empty values in Text fields?

texts:
I have written a script that checks populated (or empty) values off of coded value and range domains. It works well with the range domains. It finds empty fields as well as values that are out of range for that field. It also works well when there are values in TEXT fields that don't conform to the accepted domain codes. However, it does not find NULL (None or "") values in TEXT fields.

Here's the code.

if ftype in txtList:  ##Start coded value domain checks
with arcpy.da.SearchCursor(tbl,['UID',fname,'QAQC_FLAG']) as cur:
    listValues = []
    for row in cur:
        uid = row[0]
        fieldValue = row[1]
        if fieldValue == None or fieldValue == "":
            insRow = ["VNULL",fname,"","Value not populated with incoming data",uid]
            insCur.insertRow(insRow)
            msg = "Added row {0} to {1}.".format(insRow,qaTbl)
            arcpy.AddMessage(msg)

I can't figure out why it won't identify empty fields in the table. Some of these fields have ALL empty values and the script isn't identifying them.


texts: Android API - how get the car position (fixed to the road)?

texts:
All the android.location classes receives best REAL position of the GPS.
Is that any way to get from Google navigation API the CAR position x,y - fixed (tighten, interpolated) to the closest road on the Google map?


texts: how to calculate RMS 2D and RMS 1D, using python?

texts:
I would like to calculate:

RMS 2D: standard deviation for the two coordinates of flat (67% confidence level)
and
RMS 1D: standard deviation for height (67% confidence level)
What is the formula and what parameters do I need for dual frquency survey grade receiver?
At the moment I haves something like:

import math
spaceSegment = 0.4
controlSegment = 3.5
tropo = 0.2
iono = 0.4
multipath = 0.6
receiver_error = 0.3

hdop = float(raw_input("Enter hdop "))

vdop = float(raw_input("Enter vdop "))

ure = math.sqrt(spaceSegment**2 + controlSegment**2)
uee = math.sqrt(tropo**2 + iono**2 + multipath**2 + receiver_error**2)
uere = math.sqrt(ure**2 + uee**2)

horizontal_error = 2. * hdop * uere    #2D horizontal position error
vertical_error = 2. * vdop * uere    ##Vertical Position Error

...but I cannot find gps error budget for survey grade receivers nor formulas
for 67 and 68 percentage calculation.
Is any simple way for my computation?
Above computation is taken from matlab gps toolkit I rewrited in python but it refers to gps budget for single frequency receiver.
I also have two formulas:

and

and their decription:
Compute the rms of the radial errors, i.e. the linear distances between the measure and known (or mean) positions. It can be predicted by using covariance analysis by multiplying the HDOP, a measure of the satellite geometry, by the standard deviation of the observed pseudoranges. This predictability makes it a much more convenient measure in practice.
...so then I am completely confused and any help will be appreciated.


texts: Is it possible to change extension name TIFW to TFW in QGIS?

texts:
In QGIS from the "project-->Save as image" dropdown it is possible to save the current map canvas as Tiff. A world file is also automatically exported. Nice!
BUT, I use the TIFF+TIFW in another program that only accepts .TFW file. This is the same file, but with a 8.3-type extension.
Is it possible to get QGIS to export with the TFW extension instead of TIFW?


texts: How to Get Access to a Map been used in ArcGIS Online App

texts:
Can someone tell how I can access to GDB/ ArcGIS Online Map/ or ArcGIS Server service of the Map been used in this App (https://www.arcgis.com/apps/Cascade/index.html?appid=e8a78203d79042acbec87255d45e6ab9)?
    <iframe width="100%" height="800px" src="https://www.arcgis.com/apps/Cascade/index.html?appid=e8a78203d79042acbec87255d45e6ab9" frameborder="0" scrolling="yes"></iframe>

I need to get the Network assets/ Not sure this is been made using  Utility or Gematric Network.


texts: How do I change style/colour of a few selected features?

texts:
We have features which represent the locations of pupils within the local area and would like to change the colour of these features to represent whether they reside within the catchment area of a particular school or not?

Can anyone advise?


texts: Getting the latitude and longitude coordinates (to plot 2D) from a well survey

texts:
I basically want to do what you did here How to create a layer of directional/horizontal wells?.

I am able to get the coordinates from a well survey to create a shapefile to be able to visualize the "worm track" and from the coordinates get the middle point of the lateral length.


texts:
I open MSYS from the OSGEO4W and type wget and I get 'command not found'. I'm using Windows 7.

Does anyone know how to solve this?


texts: Setting symbology type and color based on data using QGIS

texts:
I have a point vector layer with several attributes in QGIS Dufour.
One of them is satellite name, and the other is date. I want to display different dates with different colors (using some gradient) while the satellite name rules the type of symbol.
For example GOES-13 is a square, AQUA a dot, TERRA a triangle.

Satellite  date

goes-13   01/01/2014

goes-13   01/01/2014

aqua      01/01/2014

aqua      01/01/2014

terra     02/01/2014

... and so on

I know how to change the size or rotation of the symbol, but not the type.

the solution


texts:
In QGIS2threejs,I want to change the fonts in 3D presentation.I learned that I can do it by change it in the css file through
 Different fonts in QGIS2threejs when I change the color and font,it doesn't work,the color and font doesn't change in 3D presentation.


texts: I had an excel file marking points on my map. I then added a base map and my points are gone. Even with the base map layer turned off the points do not show. I tried reordering the layers and tried deleting then re-adding the excel file and nothing seems to work. Any ideas?


texts: Python if/then Statement in Field Calculator to return 0 or 1 based on numbers divisible by 28?

texts:
I don't use python all the time.

In a point feature class, I'm trying to populate an ID field with 1's and 0's based on whether the FID is divisible by 28.

I've started the following code but can't get it to execute:

def myCalc(FID,ID):
  if FID % 28 == 0:
    return 1
  else:
    return 0


texts: Using cloud confidence to create cloud mask from Landsat 8 BQA?

texts:
I want to use the Landsat BQA (TOA images) to create a cloud mask in GEE. In particular I want to mask out all areas where there is high confidence that cloud cover is present. I have tried to adapt the example script, which utilises the single bit cloud band:

var cloudmaskL8 = function(image) {
   var qa = image.select('BQA');
   var mask = qa.bitwiseAnd(1 << 4).eq(0);
   return image.updateMask(mask);
}

But I do not understand how the 'bitwiseAnd' function works in this context and
I'm therefore unable to adapt it for the double bit confidence bands.

How does the 'bitwiseAnd' function work in this context?


texts: Joining and Calculating using ArcObjects with vb.net

texts:
I'm trying to join an excel table to a shapefile and then calculating those fields using vb.net. The code I have written has no errors, but won't work when I try it in ArcGIS 10.1 and I'm out of ideas.

Can anyone find an error in my code or direct me to other sources?

ModuleQ3 Code:

Public Sub JoinCalc(ByVal excelXLS As String, ByVal fromExcelTable As String, ByVal toFeatureClass As String, ByVal toField As String, ByVal fromField As String, ByVal calculateField As String, ByVal calculatorExpression As String)
    Dim pMxDoc As IMxDocument
    Dim pMap As IMap
    Dim pFeatureLayer As IFeatureLayer
    Dim pFeatureClass As IFeatureClass
    Dim featureIndex As Integer

    pMxDoc = My.ArcMap.Application.Document
    pMap = pMxDoc.FocusMap

    featureIndex = FindLayer(toFeatureClass)

    If featureIndex = -1 Then
        MsgBox("The join feature class is not in the current map document.")
        Exit Sub
    End If

    pFeatureLayer = pMap.Layer(featureIndex)
    pFeatureClass = pFeatureLayer.FeatureClass

    Dim pWorkspaceFactory As IWorkspaceFactory
    pWorkspaceFactory = New ExcelWorkspaceFactory

    'open excel workbook
    Dim pFworkspace As IFeatureWorkspace
    pFworkspace = pWorkspaceFactory.OpenFromFile(excelXLS, 0)

    'open excel table
    Dim pFromTable As ITable
    pFromTable = pFworkspace.OpenTable(fromExcelTable)

    'Join the table to the feature class.
    Dim pMemRelFact As IMemoryRelationshipClassFactory
    Dim pRelClass As IRelationshipClass
    Dim pDispRC As IDisplayRelationshipClass

    ' Create a memory relationship class.
    pMemRelFact = New MemoryRelationshipClassFactory
    pRelClass = pMemRelFact.Open("Table-Layer", pFromTable, fromField, pFeatureClass, toField, "FeatureClassName", "TableName", ESRI.ArcGIS.Geodatabase.esriRelCardinality.esriRelCardinalityOneToOne)

    ' Perform a join.
    pDispRC = pFeatureLayer
    pDispRC.DisplayRelationshipClass(pRelClass, ESRI.ArcGIS.Geodatabase.esriJoinType.esriLeftInnerJoin)

    MsgBox("The join is complete")

    Dim pDisplayTable As IDisplayTable
    pDisplayTable = pFeatureLayer

    Dim pFields As IFields
    pFields = pFeatureClass.Fields

    Dim fieldIndex As Integer
    fieldIndex = pFields.FindField(calculateField)

    If fieldIndex <> -1 Then
        Dim pCursor As ICursor
        pCursor = pDisplayTable.SearchDisplayTable(Nothing, True)

        Dim pCalculator As ICalculator
        pCalculator = New Calculator
        With pCalculator
            .Cursor = pCursor
            .Expression = calculatorExpression
            .Field = calculateField
        End With

        pCalculator.Calculate()

        MsgBox("It should be calculated.")

    Else
        MsgBox("The field does not exist")
        Exit Sub
    End If
End Sub

These are the parameters that I am passing through:

ModuleQ3.JoinCalc("ONT_Climate_Normals.xlsx", "ONT_Annual", "ONT_Stations_Shape", "CLIMATE_ID", "CLIMATE_ID", "TotalAnnua", "NORMAL_52 + NORMAL_54")


texts:
Does anyone have experience subsetting/clipping a large NAIP imagery file in MrSID format? Some of the NAIP files we have are by county and they have really large file sizes.

The files are from the FSA NAIP Imagery.

How does it work if you subset a MrSID file? Does ArcGIS has to uncompress it first before it can proceed to clip or subset an area of interest?


texts: Removing fields from output FeatureSet in ArcPy for geoprocessing service?

texts:
The geoprocessing service runs a script that is trying to output a FeatureLayer or a FeatureSet that will contains a selection with only some fields of the original FeatureClass.

I dont want to copy the FeatureLayer into a new FeatureClass, because I need to preserve selected features OBJECTID.

I have try to use the FieldInfo paramater to put some fields to "HIDDEN" but in the REST output those fields appear anyway.

As I only need the output in json format, my second attempt was to load the FeatureLayer in a FeatureSet, then, I remove the unwanted fields directly by formatting the json (see below).

import arcpy
import json

codeCs = arcpy.getParameter(0)
outFields = arcpy.getParameter(1)
expression = u"CODECS = '{}'".format(codeCs)
result = arcpy.MakeFeatureLayer_management(SOURCE_FC, "temp_lyr", where_clause=expression)
fSet = arcpy.FeatureSet()
fSet.load(result)
desc = arcpy.Describe(fSet)
esriJson = json.loads(desc.pjson)

fields = esriJson["fields"]
for field in fields[:]:
    if field["name"] not in outFields and field["name"] not in SOURCE_REQUIRED_FIELDS or field["name"] in SOURCE_HIDDEN_FIELDS:
        fields.remove(field)
        for idx, feat in enumerate(esriJson["features"][:]):
            del esriJson["features"][idx]["attributes"][field["name"]]

arcpy.SetParameter(2, json.dumps(esriJon))

This solution works, but it seems over-complicated to me. Anybody have an idea how I could remove some fields from a FeatureLayer?


texts: QGIS attribute form table of contents highlighting is missing

texts:
When in form view in the attribute table, QGIS displays a "table of contents" beside the form, which usually highlights the feature being displayed in the form. Example:

In two of my layers, somehow that highlight does not show. Example: (the row directly beside LongLabel should be highlighted).

Any idea what setting I've munged somewhere, and how to fix it?

The TOC is set to display a virtual field, but the problem on those layers doesn't get fixed if I use a regular field or a custom expression. I do have conditional formatting conditions (and they work), but the problem persists when I remove them too. The form is Drag and Drop rather than autogenerate, but I can't imagine drag and drop wouldn't allow highlighting, and I think it worked previously.

(Using QGIS Win 3.12, but this has been going on for some time and with several 3.x versions, I just haven't had time to try to fix it.)


texts: Convert geodatabase categorized style to QGIS style

texts:
I'm using ArcMap 10.2.2 and QGIS 2.6.1. Brighton.
I set up a geodatabase in ArcMap and used subtypes while I was digitilizing geometries.

As I'm opening the geodatabase in QGIS, I noticed that the programm is showing me all geometries in one color, and I'm not able to under categorize the styles.

Is there a solution?


texts: Visualising LiDAR data in 3D as ArcGIS Service in browser with ArcGIS Server?

texts:
I have lidar data which I would like to present over the intranet in 3D. I have arcgis server 10.3.1 with the image extention. I am however not sure if such a service is available.

Is this possible?

As far as I can see this is only available online through webscenes.


texts: How do I only view selected rasters from a geodatabase raster catalog in ArcGIS 10.0

texts:
I am using a geodatabase raster catalog to hold raster files for multiple countries, mostly for the sake of organization. I don't actually have the need (often) to visualize or draw these all at once.

When working on some country specific task, I'd like the current map document to only show the related raster. Also, for spatial analysis, I would want the selected raster to be the only raster operated on.

How do I accomplish this?


texts: Unable to create geometric network on NHDFlowline - invalid enabled field type

texts:
I need to create a geometric network in order to perform a downstream trace on a NHDFlowline data in my ArcGIS.

I have two flowline layers on one of them I am able to create a geometric network and on the other I am not. I have noticed that the difference between this two layers is in the column called "Enabled".

In the working layer this field is an short int type with "domain enabled", the not working layer has simple text enabled field (see the attachement)

I really don't know what is all about. I have no idea how can I change the field type to be able to create a geometric network.


texts: Creating boxplots from two rasters using the Raster package in R, using only a subset of classes from one

texts:
In R using the raster package, you can create box plots from two rasters (x and y) like so:

boxplot(x,y)

where:

x = Raster* object
y = If x is a RasterLayer object, y can be an additional RasterLayer to group the values of x by 'zone'

Which works perfectly when I want to group the values of x by all the classes in my thematic raster y.

But if I only want to group the values in x by a subset of the thematic classes in y, is there a way to do this?  I could replace any values in y outside the values I am interested with NODATA values, but I was hoping there a cleaner, simpler method.


texts: How to get coordinates by click on map and send coordinate to MySQL?

texts:
I'm building my own web map.

Now I want to get the coordinates by click on the map and send the coordinates to MySQL.

I tried OpenLayers to get the coordinates from Openstreetmap,but I didn't know how to send it to my MySQL.

Can anyone help me?


texts: R raster Package Moran's I interpretation

texts:
I am using raster package to calculate the local Moran's I . The example gives the range of Moran's I between - 1 to 2.47. On my own data I have see the value range -3.070423 - 7.228558  How can Moran' I value be larger than 1 ? Most of the literature point out the value of global Moran's I is between -1 and 1.
What does local Moran's I value larger than 1 mean as implemented by the R raster package?

#data
r <- raster(nrows=10, ncols=10)
r[] <- 1:ncell(r)
plot(r)

#local Moran's I
x1 <- MoranLocal(r)
plot(x1)


texts: Calculating mean value for all pixels within given distance for all pixels in study area using raster analysis in QGIS

texts:
I am working in QGIS
I have a raster layer depicting the population density for a country. A portion of this raster layer is used as a study area.
I want to know for each pixel in that study area what's the mean population density value for all points within say, 15 km of that pixel.
Where do I start?


texts: Seeking ArcGIS Pro tool to calculate network distance from point to nearest line?

texts:
I am trying to calculate the network distance from a point to the nearest poly line feature as pointed in the image below.

Is there any ArcGIS Pro tool that can calculate network distance?


texts: I'm looking for a guide on how to get up and running with ruby and gdal on Mac. I've seen some examples on how to use gdal in ruby but never found any tutorials on where to download and install it.


texts:
I want to filter a OSM file by several criterias by combining several AND and OR operators via Layer Properties -> Source -> Query Builder. My filter looks like:

"highway" in ('secondary','tertiary','unclassified','residential','secondary_link','tertiary_link','living_street','service','services','road','rest_area','yes')
and ("waterway" IS NULL or "waterway" in ('no'))
and ("barrier" IS NULL or "barrier" in ('no'))
and ("maxspeed" IS NULL or "maxspeed" in ('5','10','15','20','25','30','35','40','45','50','55','60') or "maxspeed" <= 60)
and ("access" IS NULL or "access" in ('allowed','*','public','residential','services','unknown','yes','bus'))
and ("tunnel" IS NULL or "tunnel" in ('no'))

However, it returns all features. Maxspeed currently is a String field, but thats not always the case.

How do I set up a complex filter expression for QGIS 3 correctly, or to be more specific, how do I combine several AND and OR operators?


texts:
So i"m basically getting this JSON returned from a REST call (I'm using Python).

[
  {
    "ipAddress": "X.X.X.X",
    "score": "1000",

  }
]

I want to convert this to geoJSON so I can plot the data on a leaflet map. So I plan to send this object off to freegeoip.net and get the latitude / longitude data.

But how do I convert this to geoJSON after while retaining the ip address and score data? Is there a library or something that I can use for this purpose? Any ideas would be great.


texts: API for time and distance between two addresses?

texts: The Google TOS appears to limit the use of API for time and distance between two addresses to use in conjunction with a displayed map. Does anyone know of any other mapping data source that offers a similar API without that limitation? I need JUST time and distance, for use in driver scheduling and management, not directions and such.


texts: Does uninstalling ArcGIS Desktop delete all user data as well?

texts:
Where can I get my shapefiles, layers, geodatabases etc. that were created with ArcGIS 10.1 AFTER I uninstall ArcGIS for desktop 10.1?

Are the files saved in Documents folder remaining after uninstall arcgis 10.1?


texts: Exporting Vehicle Routing Problem Solution Tables in ArcGIS ModelBuilder?

texts:
I want to export the "Orders" and "Routes" tables/shapefiles that come as a result of solving a Vehicle Routing Problem (VRP) in ModelBuilder.

However, I don't see where these layers are temporarily saved so that I can export them and store them in a specific location in my documents folder.

I have tried with "Copy Features" and "Feature Class to Feature Class" with no success (ArcGIS 10.2.2, advanced license).

Is this possible?


texts: I am trying to follow the instructions given here, Tutorial: Making Heatmaps using QGIS and GRASS ,  to build a heat map with my own data. I taken the first step and imported a csv data table into the layer bar, but when I attempt to save it as a shapefile, QGIS tells me that is not a valid or recognized data source. What am I doing wrong?


texts:
I have a footprint of some buildings (see the picture).

The selected buildings have an area below 30 sqm.
As you can in the picture see some buildings share almost borders with a building bigger than 30 sqm. When others don't.
The building below 30 sqm but shares a border with a bigger building I want to include in my calculation. The others below 30 sqm should be sorted away.

I have succeeded with this through desktop GIS, and the solution looks like this:

But I want to do the same process in postgres.

I can select the buildings pending on their area and I have put them in to two new tables:

Create TABLE building_below as
 Select * from building where area<30
and
Create TABLE building_above as
Select * from building where area>30

I can not figure out how to select layer by location in PostgreSQL.

I have tried with inspiration from:
Join polygons by location PostgreSQL/PostGIS
Merging adjacent polygons in shapefile that has been split at tile boundaries?

But it didn't work.

I have tried this:

SELECT DISTINCT ON (s.gid) s.gid,*

 FROM building_below s

LEFT JOIN building_above h ON ST_DWithin(s.geom, h.geom, 0.01)

ORDER BY s.gid, ST_Distance(s.geom, h.geom);

But it just select all the rows/buildings in building_below.

I know from desktop GIS, that it should not do that.


texts: Unable to update time frame on layer in ArcPy?

texts:
I'm working on a script to ingest some data, create some data layers with it, and publish that to ArcGIS Server.
Most of that is done with no problem, but I got a new requirement today, time-enabling some of those feature layers.  That looks pretty straight-forward, but it isn't working as I'd expect, so I'm hoping someone or some ones here can help.
# dFrame is the data frame, lyrName is the name of the new layer
def createLayer (self, dFrame, lyrName):
    lyrFileName = ...
    timeLyrFileName = ...

    arcpy.MakeFeatureLayer_management (lyrName, lyrName + self.lyrSuffix)
    arcpy.SaveToLayerFile_management (lyrName + self.lyrSuffix, lyrFileName)
    newLyr = arcpy.mapping.Layer (lyrFileName)
    ...
    # set some other layer properties that are working fine
    arcpy.mapping.AddLayer (dFrame, newLyr)

    if not os.path.isfile (timeLyrFileName):
      timeLyr = arcpy.mapping.Layer (timeLyrFileName)

    if not timeLyr is None:
      arcpy.mapping.UpdateLayerTime (dFrame, newLyr, timeLyr)
    return

(The mxd holding the data frame is saved after this method returns.)  If I put some print statements after that call to UpdateLayerTime, everything looks good.  newLyr.time.isTimeEnabled is True, newLyr.time.timeOffset is the same as timeLyr.time.timeOffset (though, oddly, both show as 0.0 years rather than 1.0 hours.  Hoping just a small bug in EsriTimeOffset's toString method (or whatever the python equivalent of that is)), and newLyr.time.timeStartField also looks correct.
The fly in the ointment now is that, when the program finishes, if I open the mxd in Desktop, the layers created are all there, but none of them are time-enabled. (FWIW, we're running the ArcPy that came with Server 10.6.)
Anyone got a clue what's going on?


texts: Adding ST_GEOMETRY column to existing table that already has X Y columns?

texts:
I have created and Oracle 11g geospatial geodatabase and have create an arcsde connection.
I have been connected using arcmap, as a simple user (not sde or schema-user), who has all the proper rights to edit tables etc..
I use arcgis for desktop 10.2.
In my geodatabase I have created all my business tables needed for my applications.
In this geodatabase there is on table tha contains X Y points.
I suppose that I need to add an st_geomery column so as to enable arcgis realise that it contains geometries (points), have indexing etc.

One option is to alter the table via an oracle query and add an st_geometry column.
In this case I cannot see the table in the listed tables in arcmap. I dont know why.
I suppose that a second option could be to edit this table through arcmap and add a st_geometry column of points. But I cannot find any option in arcmap to do this.

When I double click on my table i can see these options (missing st_geometry option)

Any help on these two options?


texts: Selection of rainfalls seasons for Sentintel-2 data in Google Earth Engine

texts:
I want to filter the Sentinel ImageCollection for the rainseasons of my ROI (March to may AND October to December) for all data available (2015-2018).
With .filterDate(startdate, enddate) I can only filter one season of one year. With my code I get an Error: Layer 1: Layer error: User memory limit exceeded. Is there another way?

var s2 = ee.ImageCollection('COPERNICUS/S2').filterDate('2016-03-01', '2016-05-31')
.filter(ee.Filter.calendarRange(3,5,'month'))
.filter(ee.Filter.calendarRange(10,12,'month'));

var composite = s2.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
              .map(maskS2clouds)
              .median();

function maskS2clouds(image) {
  var qa = image.select('QA60');
  var cloudBitMask = ee.Number(2).pow(10).int();
  var cirrusBitMask = ee.Number(2).pow(11).int();
  var mask =
  qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
  }

Map.setCenter(37.5211, 1.3182, 13);
Map.addLayer(composite, {bands: ['B3', 'B2', 'B1'], min: 0, max: 0.3});


texts: Overview of GPS tracking system technologies

texts:
Background: I've been asked by a friend to estimate the cost of setting up a GPS tracking system for a small fleet of buses, so we can do real time reporting of bus locations via a webmap.

In my research, I saw three kinds of active GPS trackers:

Radio
SMS
Satellite

Some major considerations are

minimizing long term costs, and
making the system as reliable as reasonably possible.

Secondary concerns are minimizing up front cost, latency, difficulty of installation/setup, level of support required after the initial setup is complete.

The buses run within a 5 mile radius, so radio is feasible.

With an eye towards reducing long term costs, radio seems good: it would be the only way to run the system without a middleman.

SMS requires cell phone towers, and satellites require a space program.  I am concerned that a DIY solution with radios might prove to be less reliable, though.

But if we bought SMS transmitters that provided the data through a third party API, they could go out of business, or have unacceptable down time or latency.

My question(s):

What would each of these three solutions look like?
What hardware is involved?
What are the trade-offs between them?
What problems have you encountered?
Have you used a system that Just Works?


texts:
I'm brand new to GIS (and this site) so bear with me if there is an obvious answer for this but I've done a lot of looking and haven't found anything that worked.

I have two polygon layers that show the same waterway buffer for a large area but one is more extensive than the other. I want to subtract the smaller layer from the larger one so that I am left with the areas that are only covered by the larger. This seems like the type of problem that could be easily solved with the "erase tool" but I don't have access to that. The "select by location" seems to be a popular suggestion but I've played around with it and haven't gotten it to work. Also the two layers are not quite identical.

What should I do?


texts: Placing the label style alongside the layer style within the legend - Arcmap 10

texts:
I've been searching for a way to include the label symbology within the legend, something similar to the image below: (Apologies for the quality)

There was a post some 4 years back that touched on the subject within Arcmap 9. Hoping someone has found a way to create this since then.

Anyone have a workaround to create something like this? Other than creating a text box for each layer?


texts: Getting WMS server from GLOBELAND30 Landcover interactive map?

texts:
I am very interested to get GLOBELAND30 WMS server. To get it I've used some hacking:

1) I've inspected it with Chrome Inspector and did some actions (see below):

2) Got WMS server (see below)

3) Added the WMS url http://218.244.250.80:8080/erdas-apollo/coverage/CGLC?LAYERS=cglc30_2010_0 to QGIS where it shows WMS layers but can't load the data.

Is there any idea how to solve this?


texts: How to select points within a certain distance from a line?

texts:
I have a file with a points layer and a lines layer. I would like to select only the points that lay within certain distance (a couple of meters) of the lines.

Is there an easy way to do this?

Or should I create a new polygone layer wich forms a kind of buffer around te lines and then use the select bij location tool?

How can I create such a layer?


texts: JTS - Java Topology Suite - Calculate radius from a given point

texts:
From point A, I want to check if a point B is within a certain radius from A.
How calculate it usindo JTS?

I tried this:

Coordinate posicaoAtual = new Coordinate(-12.974900375893823, -38.48034163207694);
Localidade localidade = new Localidade(-12.976353604547635, -38.47492356986686);
Geometry ponto = new GeometryFactory().createPoint(localidade.getCoordenada());
Geometry pos = new GeometryFactory().createPoint(posicaoAtual);
Geometry raio = ponto.buffer(100.0, 8); // HERE
System.out.println(raio.contains(pos));

My doubt is to know which the distance unit from the first parameter of the buffer method?


texts:
I am using ArcMap 10.5.
In the green polygon of the area below there are some small gaps between the hatches. It is interesting that these gaps appear in the pdf exported from the ArcGIS (we used 450 dpi) and in the layout view as well if you use a similar scale as the one showed in the image.
If you zoom in, these gaps between the hatches disappear.
Does someone know a solution for this?
I took a photo of the printed version of the map:


texts: Setting projection parameters in Geopandas

texts:
I have a shapefile which I want to reproject to lambert conformal conic and set the projection parameters such as the following:

False_Easting, False_Northing, Central_Meridian, Standard_Parallel_1, Standard_Parallel_2, Latitude_Of_Origin

I used to do this easily with ArcMap, but I do not have any access to ArcMap anymore. Also I can't seem to do this with reproject under the vector toolbox in QGIS.

How could I do this in geopandas?


texts: How to correctly resize raster GIS images to a given px width?

texts:
I work on ETOPO1, a 466MB file, and want to move to SRTM (~20GB Earth). My workflow involve several time-consuming processing steps so I strongly want to work on a minimum raster. It usually want to take 10.000 to 100.000px width GIS rasters to scale them down to ~1000px width, like my final map.

How to correctly resize a raster GIS file to a given width via terminal ? (gdal prefered)

It should keep negative and positive z-values with precision.

Note: I already tried the following...

Crop: to crop my area of insterest using gdal_translate -projwin $(WEST) $(NORTH) $(EAST) $(SOUTH) ETOPO1_Ice_g_geotiff.tif crop.tmp.tif, it still give me files 50 times too big for a processing 50 times slower.
Early resizing: early resizing with Imagemagick + reinjecting the GIS metadata for following processings destroys the negative values (altitudes under sea level become 0m).


texts: How export geodatabase version differences in ArcPy?

texts:
With ArcObjects I have this code in C# to export the differences to a XML file:

        ESRI.ArcGIS.GeoDatabaseDistributed.IVersionDataChangesInit pVDataChgs;
        var pvTemp = new ESRI.ArcGIS.GeoDatabaseDistributed.VersionDataChanges();
        pVDataChgs = pvTemp as ESRI.ArcGIS.GeoDatabaseDistributed.IVersionDataChangesInit;
        pVDataChgs.Init(pwkspname, pParWkSpname);
        var pvDataChanges = pVDataChgs as ESRI.ArcGIS.GeoDatabaseDistributed.IDataChanges;

        //Export to an XML file
        ESRI.ArcGIS.GeoDatabaseDistributed.IExportDataChanges pDCExp;
        pDCExp = new ESRI.ArcGIS.GeoDatabaseDistributed.DataChangesExporter();
        pDCExp.ExportDataChanges(sXMLExpFile, ESRI.ArcGIS.GeoDatabaseDistributed.esriExportDataChangesOption.esriExportToXML, pvDataChanges, true);

I need an equivalent functions to get the same result trought ArcPy. But the closest thing I've found is this function: FeatureCompare_management.

How can I compare versions, not features one by one?


texts: GeoTIFFs made in GDAL don't have min/max range when loaded in Arc/QGIS

texts:
I'm creating some floating point rasters in GDAL with the GeoTIFF driver.  When I load the resulting images into QGIS or Arc the default symbology puts the min at -3.40282e38 and max to 3.40282e38, so the raster looks exactly gray.  Is there a way I can write the actual range into the GeoTIFF directly so when its loaded into a GIS program it automatically scales the histogram nicely?  I've tried creating a default histogram like this:

rasterMin, rasterMax = raster.GetRasterBand(1).ComputeRasterMinMax()
raster.GetRasterBand(1).SetDefaultHistogram(rasterMin, rasterMax, 255)

But I get a error that looks like a binding issue:

  File "/usr/lib/python2.6/dist-packages/osgeo/gdal.py", line 846, in SetDefaultHistogram
    return _gdal.Band_SetDefaultHistogram(self, *args)
TypeError: not a sequence

Any suggestions on what I could do, or if I've made a mistake?


texts:
Is it possible to return null voxels (those without any points) using voxel_metrics()?

As the description is written, voxels are only returned when they contain at least one point. I am currently attempting to look at ratios of empty and not empty voxels in various height layers such as discussed in Pearse et al. (2019).

Reference:

  Pearse, G. D., Watt, M. S., Dash, J. P., Stone, C., & Caccamo, G.
  (2019). Comparison of models describing forest inventory attributes
  using standard and voxel-based lidar predictors across a range of
  pulse densities. International Journal of Applied Earth Observation
  and Geoinformation, 78, 341-351.


texts: Integrating raster datasets using ArcGIS for Desktop?

texts:
I have two raster datasets. One raster dataset shows population density, with the range of population density values given a score from 1 to 10. The score is entered into a separate column.

The second raster dataset is average income. The range of average income values have been given a score from 1 to 10. The score is entered into a separate column.

I would like to integrate these two raster datasets and create a third. I would like the third raster dataset to display the sum of the first two raster dataset scores. So the third raster dataset will have values ranging from 2 to 20.

I am using ArcMap 10.3 with the spatial analyst extension.

Edit: The context is conservation planning. We have a number of layers that we would like to integrate to develop an overall "cost score". The individual input layers have been scored so that a value of 10 indicates highly modified or impacted / not worth protecting. A lower score indicates the opposite. This will feed into a larger analysis involving species' distribution, endemicity, etc.


texts: indexes to optimize ST_Contains and ST_Overlaps for geoms in two tables

texts:
I'm a programmer getting into the wonderful world of GIS. I just had a fun day optimizing queries (took one from 60s to 3.4s, and another from 38s to 1.2s). Below is the last query that I can't seem to optimize! (My white whale!)

my query is basically the following

SELECT child,st_asgeojson("child".wkb_geometry)
FROM child
WHERE ST_Contains(
(SELECT wkb_geometry FROM parent WHERE ogc_fid=xxx),
"child".wkb_geometry
) OR ST_Overlaps(
(SELECT wkb_geometry FROM parent WHERE ogc_fid=xxx),
"child".wkb_geometry
);

I understand that this OR is probably terrible, but truth is the query is terrible with and without it... (both pieces are individually terrible). It's currently at ~40s average. I'm fine with optimizations that work for one or the other piece of the query...(or both separately, or both together).

I just learned about indexes and the last 2 queries I optimized got better with basic btree indexes (and one multi-column index)...I assume I need a combination of an expression index and something else, but I don't know where to start thinking about the solution. I also don't know how to have a cross-table index or if that is necessary...for this query.


texts:
I have created my own script in ArcGIS 10.1 and for one of my parameters I have a drop down list that the user has to select from.  The longest entry in that list is only 6 characters, but the width of the drop down box is the total width of the window.

Is there a way to set the drop down box width to something smaller?


texts:
I have a polyline drawing from AutoCAD that I would like to reposition manually, if possible. During my efforts so far it has been turned into both a tab file and shapefile.

The main source of difficulty appears to arise from the fact that it was imported without a CRS and when I force one onto it it disappears.

Would anyone be able to render assistance?


texts: Calculate distance in meters using results in degrees

texts:
I receive some distance differences in degrees, like this:

0.0009698680322606666

using this function:

def radius(gpd_df, cpt, radius):
    sindex = gpd_df.sindex

    bbox = (cpt.geometry.x-radius,
            cpt.geometry.y-radius,
            cpt.geometry.x+radius,
            cpt.geometry.y+radius)

    good = []
    for n in sindex.intersection(bbox):
        dist = cpt.geometry.distance(gpd_df['geometry'][n])
        if dist<radius:
            good.append((dist,n))

    good.sort()

    return [x for x in good]

How do I transform that distance into meters in Python?


texts: Unable to read in shapefile using GeoPandas

texts:
I am trying to read in a shapefile:
sph = gpd.read_file("taxi_zones.shp", encoding="utf-8")

but I am getting this error:
'taxi_zones.shp' not recognized as a supported file format.

Not entirely sure what this means.


texts:
i'd like to take high quality photos with my tablet/smartphone during a survey with QFIELD.
But i have some issues because with the camera used by the App, i cannot zoom in/out and the quality is bad, while using the native android camera of the smartphone QFIELD crashes.

Any suggestion?


texts: How to remove temporary layers so they don't show up in parameter lists

texts: ArcGis 10.2, arcpy. I have several tools in a Python toolbox that have a feature class parameter.  The parameter is defined as parameterType = “GPFeatureLayer”.   When running several tools in succession the list of choices for the parameter becomes cluttered with temporary layers created from previously run tools.   By temporary I mean layers created with arcpy.MakeFeatureLayer_management().  This list will sometimes include multiple occurrences of the same layer name.  I would like to remove all layers from wherever they are coming from so only those layers that are in the TOC show up as choices in the tool parameters.  I have tried arcpy.Delete_management() but that resulted in deleting the actual features the layer was generated from rather than the layer (OUCH!).   I have also tried emptying out the scratch GDB.  Thanks in advance for your input.


texts: Licensing ArcGIS Enterprise when scaling in cloud environment

texts:
Esri's ArcGIS Enterprise licenses are for use on machines with "up to 4 cores".  I understand that in cloud environments (such as AWS) that available machines have an equivalent number of physical cores associated with them (breakdown here https://aws.amazon.com/ec2/physicalcores/).
How would dynamic scaling of computing resources work in a scenario where resources need to be scaled beyond the 4 cores that are licensed?
For example, would two ArcGIS Enterprise licenses need to be purchased if an organization needs to intermittently scale up to the equivalent of an 8 core machine in an AWS environment?


texts:
I want to remove the land area based on mask file; which is a binary image. When I use following gdal_calc.py command:

gdal_calc.py -A G:\2016.06.MOD\MODIS_JPG\Resize\Arctic.2016059.aqua.1km_Resize.tif -B G:\2016.06.MOD\MODIS_JPG\Mask_Land\Mask.tif --outfile=G:\2016.06.MOD\MODIS_JPG\Mask_Land\test.tif --calc='A*B ' --NoDataValue=0

it was produced following ImportError message:

Traceback (most recent call last):
  File "D:\PROGRA~1\QGIS2~1.14\bin\gdal_calc.py", line 50, in <module>
    import numpy
  File "D:\PROGRA~1\QGIS2~1.14\apps\Python27\lib\site-packages\numpy\__init__.py", line 142, in <module>
    from . import add_newdocs
  File "D:\PROGRA~1\QGIS2~1.14\apps\Python27\lib\site-packages\numpy\add_newdocs.py", line 13, in <module>
    from numpy.lib import add_newdoc
  File "D:\PROGRA~1\QGIS2~1.14\apps\Python27\lib\site-packages\numpy\lib\__init__.py", line 8, in <module>
    from .type_check import *
  File "D:\PROGRA~1\QGIS2~1.14\apps\Python27\lib\site-packages\numpy\lib\type_check.py", line 11, in <module>
    import numpy.core.numeric as _nx
  File "D:\PROGRA~1\QGIS2~1.14\apps\Python27\lib\site-packages\numpy\core\__init__.py", line 24, in <module>
    raise ImportError(msg)
ImportError:
Importing the multiarray numpy extension module failed.  Most
likely you are trying to import a failed build of numpy.
If you're working with a numpy git repo, try `git clean -xdf` (removes all
files not under version control).  Otherwise reinstall numpy.

Can you give me some advice for the error?


texts: Changing Default Directory for ArcCatalog Connect To Folder?

texts:
When connecting to a folder in Arccatalog, is it possible to change the directory that the Connect To Folder dialogue box defaults to?

It currently defaults to C:\Users\username


texts:
I am not even sure how to ask this question and I have Googled what I think is wrong with no success.

When I create a map with different shadings on it, it looks perfect until it is wrote to PDF.

I want the hatching/shading to display as it does on screen.

what is printed


texts:
I have downloaded the Hydro 1k dataset for Africa from the USGS Earth Explorer portal.

Rasters are provided as .bil files, and according to the Readme are in Lambert Azimuthal Equal Area projection.

I'm using QGIS (ver. 2.14.3), I would like to reproject these files to EPSG 4326 - WGS84. I can run Raster > Preojections > Warp, setting the Target SRS to WGS84 and it runs without apparent errors. However, when I load the resulting .tif into another project with some .shp files that are already in WGS84, they are at completely different scales. The Hydro 1k layers are huge, at a far large scale than the shapefiles.

Can anyone suggest what's going wrong here? Am I making any obvious mistakes? Are there settings I should alter in the Warp tool apart from the Target SRS?

EDIT

Following up on @IanTurton's comments, there is a projection file with some information I think is relevant. I'm not sure how to use this information correctly though;

PROJCS["US National Atlas Equal Area",GEOGCS["Unspecified datum based upon the Clarke 1866 Authalic Sphere",

DATUM["D_Sphere_Clarke_1866_Authalic",

SPHEROID["Clarke_1866_Authalic_Sphere",6370997,0]],

PRIMEM["Greenwich",0],

UNIT["Degree",0.017453292519943295]],

PROJECTION["Lambert_Azimuthal_Equal_Area"],

PARAMETER["latitude_of_origin",5],

PARAMETER["central_meridian",20],

PARAMETER["false_easting",0],

PARAMETER["false_northing",0],

UNIT["Meter",1]]


texts: Converting symbols from ArcGIS Desktop to QGIS?

texts:
How do I export map information display styles from ArcGIS Desktop to QGIS?

I tried exporting MXD to QGS (using mxd2qgs) - but the styles were not saved;Also, when I try to export styles from arcgis to SLD format (using ArcGIS_SLD_Converter) and finally got:

  "Index outside the range."

The index should be a positive number, and its size should not exceed the size of the collection
Parameter name: index
in System.collections.ArrayList.get.item (lnt32 index)
in ArcGis_SLD_Converter.Output_SLD.WriteToSLD()"


texts: Using a Wacom Tablet and pen for QGIS, can I pause my free-hand polygon before I am finished, and then continue?

texts:
Is there a way to pause the free-hand tool in QGIS so that I can lift the pen off of the screen, and then continue my polygon again when I am ready?
I have learned how to pan while free-hand drawing but I'd like to be able to take my pen off the screen, maybe have it create a node point without closing the polygon, and then being able to go back and continue from that node to complete my polygon?
Does such a thing exist?


texts:
I'm in the process of converting an older ArcMap extension to an Add-In in C# and I am trying to best take advantage of the XML based configuration.

Config.esriaddinx works great for the basic set up (toolbar, buttons, extension info, etc), but is there a way to add other configuration parameters, such as field names, layer names, database connections, etc?

Would I need to create an App.config file?  If so, would the App.config file still be packaged/deployed the same as the rest of the add-in?  What would be the best method to read the parameters?  I am a newbie to add-ins and any suggestions or samples would be greatly appreciated.  Thanks!


texts: Can't get the coordinates from OpenLayers polygon drawing

texts:
I'm building an app using Ionic 4 included OpenLayers.
After drawing the polygon, I can see the coordinates on the console.log but I don't have a solution to get these coordinates.

Here is what I tried

ngOnInit(){
    const scope = this;
   this.addInteractions(GeometryType.POLYGON);
   this.source.on('addfeature', function(coord){
    var feature = coord.feature;
   var coords = feature.getGeometry().getType;
  }

    ----

 addInteractions(geoTy) {
  this.map.getInteractions().clear;

  this.draw = new Draw({
  source: this.source,
  type: geoTy
});

this.map.addInteraction(this.draw);
this.snap = new Snap({source: this.source});
this.map.addInteraction(this.snap);

console.log(this.map.getInteractions().get);
 }

I haven't added things that are less important for this.

Ionic:

Ionic CLI : 5.4.14 (C:\Users\FQDev\AppData\Roaming\npm\node_modules\ionic)
Ionic Framework : @ionic/angular 4.11.8
@angular-devkit/build-angular : 0.801.3
@angular-devkit/schematics : 8.1.3
@angular/cli : 8.1.3
@ionic/angular-toolkit : 2.1.2

Cordova:

Cordova CLI : 9.0.0 (cordova-lib@9.0.1)
Cordova Platforms : not available
Cordova Plugins : cordova-plugin-ionic-keyboard 2.2.0, cordova-plugin-ionic-webview 4.1.3, (and
7 other plugins)

Utility:

cordova-res : 0.8.1
native-run : 0.3.0

System:

Android SDK Tools : 26.1.1 (C:\Users\FQDev\AppData\Local\Android\Sdk)
NodeJS : v10.15.3 (D:\Program Files\nodejs\node.exe)
npm : 6.4.1
OS : Windows 10


texts: Algorithm to merge adjacent raster datasets

texts:
I am trying to figure out how to merge adjacent raster datasets.
Is there any clear algorithms (may be a pseudo code)?